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ABSTRACT 

An analytical technique is developed to solve nonlinear longitudinal 
combustion instability problems associated with liquid propellant rocket 
motors. The analysis produces the transient and limit-cycle behavior of 
unstable motors and the threshold amplitude required to trigger a linearly 
stable motor into unstable operation. The limit cycle waveforms are 
found to exhibit shock wave characteristics for most unstable engine operating 
conditions . A method of correlating the analytical solutions with experi- 
mental data is developed. Calculated results indicate that a second order 
solution adequately describes the behavior of combustion instability 
oscillations over a broad range of engine operating conditions, but that 
higher order effects must be accounted for in order to investigate engine 
triggering. 
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SUMMARY 


This report describes an analytical technique for the analysis of 
nonlinear longitudinal combustion instability. The applications of this 
technique, which is based on the Method of Weighted Residuals, is demonstrated 
for a liquid propellant rocket combustor with a high impedance injector and a 
short nozzle. Crocco’s time lag theory is used to describe the unsteady 
combustion response. 

The methodology developed in this report can predict the linear stability 
limits of the engine, triggering limits, and both the transient and final 
periodic behavior of the combustion instability oscillations. It is shown 
that the final periodic behavior of the instability is only dependent upon the 
engine operating conditions (i.e., the Mach number, the combustion parameters n 
and t, etc.) , and is independent of the characteristics of the initial 
disturbance. Computed results show that in most cases the final waveforms 
exhibit shock-wave-type behavior, and that the number of shocks present within 
the combustor is determined by engine operating conditions. It is also found 
that the waveform of the resulting instability depends upon the proximity of 
the operating point to engine resonant conditions (i.e., to the minimum point 
on the linear stability limit) . The predicted waveforms are used to develop 
a technique for correlating the analytical solutions with experimental data. 



INTRODUCTION 


Experimental observations of rocket motors experiencing high frequency 
longitudinal combustion instability show that in a majority of cases the 
oscillations of the gas inside the combustor exhibit shock-wave characteristics. 
These flow oscillations can be initiated in two ways. In an intrinsically 
unstable engine the instabilities are spontaneous in nature and they result 
from any perturbation of the steady state flow field. On the other hand, some 
engines require the introduction of a finite amplitude disturbance to produce 
■unstable combustion. In either case, the oscillations experience a transient 
phase prior to the establishment of stable, periodic (i.e., stable limit cycle) 
waves with oscillation frequencies that are usually close to the frequency of 
one of the chamber’s acoustic modes. These observations suggest that a 
nonlinear analysis , capable of predicting the limit cycle waveforms and the 
conditions for which unstable combustion can be initiated by finite amplitude 
disturbances of the steady state flow, is required. 

In this report, the Galerkin method, that is a special application of 

the Method of Weighted Residuals (see Ref. 1 for discussion of this method) , 

is used to develop an approximate mathematical technique for analyzing the 

nonlinear behavior of rocket engines susceptible to longitudinal mode 

combustion instabilities. The desired mathematical techniques are developed 

by investigating longitudinal combustion instabilities in liquid propellant 

rocket combustors with a high impedance injector and a short nozzle. The 

Mach number of the combustor mean flow is assumed to be small. Crocco's 

2 

pressure sensitive time lag model is used to describe the unsteady combustion 
process . 

The problem is analyzed by solving the conservation equations describing 
the behavior of large amplitude combustion-driven oscillations in low Mach 
number mean flows. Because the solution of these equations requires a 
relatively large amount of computations, a simpler but more restrictive 
second order analysis is developed concurrently. In this second order analysis, 
the amplitude of the flow oscillations are restricted to be of the order of 
the steady state flow Mach number, and terms of order higher than second order 
are neglected. Hereafter, the former problem formulation will be referred to 


2 



as the third order theory, and the latter analysis will be called the 
second order theory. The applicability of the second order theory will be 
determined by comparing its results with third order solutions . It will be 
shown that from a practical point of view, the results predicted by the second 
order theory are comparable to those found by the solutions of the third order 
theory. 

The results obtained in this investigation are used to develop a 
technique for correlating the analytical results with experimental data. An 
empirical method for predicting the nonlinear waveforms is also discussed. 

A User's Manual for the required computer programs is included in the 
appendices of this report. 


SYMBOLS 


C k (t), D k (t) 


time -dependent coefficients in the series defined in 
Eqs . (10) through ( 13 ) 

semi-empirical peak amplitude of the first harmonic 
defined in fiq. ( 28 ) 



boundary conditions defined in Eqs. ( 5 ) , (6) , (8) , and (9) 


c sonic velocity 

E l5 Eg, E , E^ flow equations, defined by Eqs. (l) , (2), ( 3 ), and ( 7 ) 

I (k,£), Ig(k,m,t) , space integrals defined in Eqs. (22) through (26) 

I (k,m,-t) , I^(k,m,t) , 

Icj (k,m,t) 


k,m,t 

L 

N 

n 

P 

q. 


summation indices and axial mode numbers 
combustor length 

number of terms used in series expansions 
interaction index 

*, * *2 

dimensionless pressure, yp /Pq c q 

dimensionless acoustic-type velocity, defined in Eq. (l4) 
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0 
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/ 
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6n 

e 
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Subscripts 

e 

k,m,b 

LS 

R 

t 

z 


dimensionless wave period, T Cq/ L 

•X- -x- 

dimensionless time, t c^/L 
time correlation parameter 

-X- -X- 

dimensionless velocity, u /cq 

unsteady combustion mass source 

dimensionless axial coordinate, z /L 

dimensionless axial coordinate used in experi m ental 
correlation 

power exponent in experimental correlation 

specific heat ratio 

dimensionless peak-to-peak amplitude 

vertical displacement at constant t in the n-f plane 

ordering parameter 

-X- 

dimensionless specific volume, v p Q 

•X - 

dimensionless density, p/p Q 

■X" -X* 

dimensionless pressure sensitive time lag, f c^/L 
velocity potential 

~X- *X- 

dimensionless angular frequency, co b/c Q 
correlation parameter 


evaluated at the nozzle entrance 
axial mode numbers 

evaluated at the linear stability limit 
evaluated at retarded time, t = t-f 
time derivative 
space derivative 
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0 


injector face stagnation quantity 


Superscripts 

perturbation quantity 
~ steady state quantity 

* dimensional quantity 

~ approximate solution 


DEVELOPMENT AND SOLUTION OF THE EQUATIONS 


Problem Formulation 

An analytical technique for investigating the nonlinear stability of 
combustion-driven axial mode oscillations in liquid propellant rocket 
combustors is developed. The combustor geometry is shown in Fig. (1) . The 
liquid propellants are injected uniformly through a high impedance injector, 
converted by a complex combustion process into hot gases, and the gas products 
are exhausted through a short nozzle. The nondimens ional coordinate system 
is defined with the origin at the injector face and the nozzle entrance plane 
at z = z /L = 1. The thermodynamic variables are normalized by the correspond- 
ing injector face stagnation quantities, the velocity is nondimens ionali zed 
by the injector face steady state stagnation sonic velocity, and time is 
normalized by a characteristic time defined as the ratio of the combustor 
length to the injector face stagnation sonic velocity. 

In order to develop a problem formulation that is both simple and 
physically meaningful, the following assumptions are made: 

1. The flow is one -dimens ional, with the velocity vector parallel to the 
combustor axis . 

2. The mean flow Mach number and its derivative are small. 

3 . The flow consists of a single constituent perfect gas and liquid 
droplets of negligible volume. 
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Distributed 

Combustion 



Injector Plate Boundary Condition: 

u'(0,t) = 0 

Short Nozzle Boundary Condition 


= constant 

c + c 


Figure 1. Combustor Geometry and Boundary Conditions 
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4 . Viscosity, diffusion, and heat conduction are negligible. 

5. The liquid droplet specific stagnation enthalpy remains constant as 
the droplets move through the combustor. 

6. The momentum sources arising from gas-droplet interaction are 

negligible . 

Even with these restrictions , the equations describing the behavior of 
axial combustion instability oscillations are quite complex . A simplification 
of the analysis results when the relative importance of the various terms in 
the conservation equations is established by using order of magnitude 
arguments. In this analysis, the magnitude of each term appearing in the 
conservation equations is evaluated and all terms whose magnitudes are 
smaller than a certain threshold value are neglected. To accomplish this 
task, two ordering parameters are used. One parameter, u g , is a measure of 
the effect of the presence of mean flow upon the oscillations. The second 
parameter, e, is a measure of the amplitude of the flow oscillations. This 
investigation is concerned with the behavior of moderate and large amplitude 
instabilities in combustors having low Mach number mean flows. Consequently, 
terms of order higher than 0(u ) and 0(u e ) are neglected. Under these 
restrictions, the equations describing the behavior of the combustion 

O 

instability oscillations can be written as follows-*: 

1. Continuity: 


E 


.1 


/ . — / 

Vj + uv 
t z 



/ / 

V u 

z 


/ / _ / / 

u + w + 2v w 
z z L z 




= 0 


( 1 ) 


2 . Momentum : 


E 2 = U t + SU Z 


, du / , / / , 1 

+ — u +uu + — 
dz z y 


/ / . 1 / n 

v p + — p = 0 
z y z 


( 2 ) 


Energy: 


E 3 = P t ' + up ^ + up' + yu' 


/ du / , / / , 

Yw + v — p +ypu + 
T z 1 dz ^ z 


v(v-i) 

- 2 


du 

dz 


u 


' 2 J = 0 (3) 
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In writing these equations, the specific volume, v, is used instead of the 
density, p, in order to simplify the numerical solution of the ordinary 
differential equations that result from the application of the Galerkin 


method" 5 . 


3w / 


The term ^ represents the unsteady mass generation "by the burning fuel. 

In the present analysis this unsteady mass source (or combustion response 

function, as it is sometimes referred to) is described by Crocco’s pressure 

2 dw ^ 

sensitive time-lag hypothesis . Accordingly, ^ is then given by the follow- 
ing relation: 


|f = n ff [p'(z>t) -p'(z,t-T)] (4) 

where n and f are the two parameters that Crocco used to describe the unsteady 
combustion process. The interaction index, n, is a measure of the sensitivity 
of the combustion process to flow oscillations and the sensitive time lag, t, 
is representative of the time required for the unsteady combustion process to 
respond to flow perturbations. 

The computed solutions must satisfy the solid wall boundary condition at 
the injector face; that is 


B 1 = u'(0,t) = 0 


(5) 


and the quasi-steady short nozzle boundary condition 


3,4 


at the nozzle entrance 


B 2 = u'(l,t) - ^ u e p'(l,t) + [u e p /2 (l,t)] = 0 (6) 

Additional simplification of the conservation equations is possible when 

the amplitude of the instability is moderate. Under this condition, it is 

possible to assume that the ordering parameters u g and e are of the same order 

of magnitude and all terms of order higher than second (e.g., terms of 

0 (u g e) or 0 (e^)) are negligible; all such terms are bracketed in Eqs, (l) 

through (3). When these terms are neglected, Eqs. ( 1 ) through (3) can be 

5 6 

combined^ 5 into the following nonlinear wave equation: 
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\ = V - '•’tt - ^zt - Y S - 2 Vzt ■ (Y - 1) % S - If = 


where 9 is the velocity potential defined hy 9„ = u ( z , t) . The bracketed term 

z 

in the short nozzle boundary condition is also neglected, and the boundary 
conditions are written in terms of the velocity potential 9 : 

B_ = 9,(0,!) = 0 (8) 

J> Z 

B 4 = 9 z (l,t) + u e 9 t (l,t) = 0 (9) 

It has been shown^ that solutions of Eqs. (7) through (9) adequately describe 
the behavior of moderate amplitude axial instabilities, but that solutions of 
Eqs. (1) through (6) are required to investigate axial mode triggering. 
Solutions of both formulations of the problem are developed in this report. 

Solution Technique 

Closed-form mathematical solutions of the equations developed in the 
preceding section are not known. As a result, it is necessary to resort to 
the use of either numerical solution techniques or approximate analytical 
techniques. The former approach is generally quite complex and its application 
requires excessive computation time; furthermore, the use of numerical- 
solution techniques in general provides little physical insight into the 
problem. An appealing approximate analytical technique has been developed by 
Zinn and Powell^ who investigated nonlinear transverse combustion insta- 
bility problems. In these investigations, the undetermined function form of 
the Galerkin method, that is a special application of the Method of Weighted 

~| O 

Residuals 5 (MWR) , is used to find the desired solutions. 

In order to use the Galerkin method, it is necessary to represent the 

dependent variables by means of approximate series expansions. The proper 

choice of the series expansion is critical to the usefulness of the Galerkin 

method. Various guide lines for the choice of the approximate series expansion 

1 8 

are offered in the literature ’ . In studies of combustion instabilities it 
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is convenient to use available experimental data, which indicates that the 
behavior of high frequency combustion instability oscillations resembles the 
behavior of the chamber’s acoustic modes, as a guide. This information 
suggests that the dependent variables of the problem be expanded in terms of 
the natural acoustic modes of the chamber; each having an u nkn own time- 
dependent coefficient. Using available acoustic solutions as a guide, the 
following approximate series representations for the dependent variables are 
used: 


N 

v'(z,t) = n(t) 
k=i K 

cos (krrz) 

(10) 

N 

p ( z , t) = £ B (t) 

cos (krrz) 

(11) 

k=l 


N 

q'(z,t) = £ C, (t) 

sin(krrz) 

(12) 

k=l 


~ N 

9 (z,t) = £ n (t) 

1 T K 

cos(kTTz) 

( 13 ) 


k=l 


the variable q' represents the "acoustic portion" of the velocity perturbation 
u that is given by 


u 


'(z,t) = [-^ u e p "(zjt) 


y 2 -l 

• — V 


8y c 


/2 (z,t)]z +q / (z,t) 


( 14 ) 


The particular choice of an expression for u 7 , as given in Eq. (l 4 ) , was 
dictated by the requirement that the dependent variable satisfy the problem’s 
boundary conditions 3 (i.e., Eqs . (5) and (6)). 

The unknown time-dependent mode- amplitudes (e.g., A (t) , C (t) , etc.) 
are determined by the following mathematical procedure. The assumed series 
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expansions are substituted into the conservation equations and boundary 
conditions to form differential equation residuals and boundary residuals. 

If a residual is identically zero, then the corresponding equation or boundary 
condition is identically satisfied. On the other hand, when the equation or 
boundary residuals are not identically zero, the residuals are the errors 
that resulted from using the approximating expansions of the dependent 
variables. According to the Galerkin method, these errors (i.e., residuals) 
can be minimized in some average sense by requiring that the residuals satisfy 
certain orthogonality conditions 5 ^ In the solution of the problem defined 
by Eqs. (l) through (3), the boundary conditions, Eqs. (5) and (6) are 
identically satisfied by the chosen series expansions, Eqs. (10) through (12), 
and the required orthogonality conditions are defined by the following 

Q 

relations : 


J E^ cos(tnz)dz =0 & = 1,...,N 


1 

J E^ sin(-tnz) dz = 0 <t = 1, . . . ,N 


[ E* cos(£rrz)dz =0 -t = 1,...,N 

J 0 3 


(15) 


( 16 ) 


(17) 


On the other hand, the expansion of the velocity potential used in the 

second order solution (i.e., Eq. (13)) does not satisfy the quasi-steady short 

nozzle boundary condition, Eq. (9) • In this case, the required orthogonality 
3 9 

condition is^’ : 

1 

J E cos ('t-TTz) dz - u e 9 t (ljt) cos (£rr) =0 l = 1,...,N (18) 

0 3 

The last term in the above equation represents the effect of the nozzle boundary 
condition residual. 
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Performing the operations indicated in Eqs. (15) through (17) yields the 
following system of quasi-linear ordinary differential equations describing 
the behavior of the unknown mode- amplitude functions: 


dA 


l 


dt~ " \ " ( * T > C t + \ B l - 5 e A t ■ nu e<V\> 


N r 

+ 2 S - Q 1 (kTT)l 1 (k,-t)B 1 


k=l 


+ 2 [(mir)l 2 (k 3 m 3 ^)C k A m + (mit) I (k 3 m 3 ^) A^C^ 


m=l 


- u e I 3 (k 3 m,^)A k A m + Q 1 (rnn)l 4 (k 3 m 3 ^)B k A m 

- Q 1 (nrrr)l 4 (k 3 m 3 ^)A k B m + Q^I (k 3 m,*) A^ 

~ 2nu e I 3 (k 3 m , 1 ) (B k -B^)A m 
+ Q 2 I 3 (k 3 m 3 ^)B k B m - 2Q 2 (mJT)l 4 (k 3 m 3 ^)B k B m ]} 


(19) 


dB^ 

w - \ = - ^ C i ~ v*e B i - + n( B r B ^) 


n r_ 

+ 2 2 -Ju e (krr)i 1 (k 3 'l)B k + ^(knJl^kj^Bj 


k=l 


+ 2 ["(ran) I (k 3 m 3 £) C,B - y(inrr)l (k 3 m 3 -t)B,C 

m=l d k m 3 k 1 


m 


- Y ^2 ^ u e I 2^ k ’ m, ^ C k C m " BjB m 

+ (Y+l)Q 1 (mTT)l 4 (k 3 m 3 -t)B k B m - yQgl (k,m 3 -t) B^ 
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+ 2yQ 2 (hut) 1^ (k ,m , l) B^B J} 


( 20 ) 


where 


dC 

dt 


- = -22 {q I (k,£)F + 2Q S I. (k,£,m)B k F 2 } + F. 

k=l k m=l m 


>1 


( 21 ) 


F = — B p - u C« + 2 £ {- u(kTT)l ,(£,k)C 

3 y l e t L e x 


N r-i 

+ £ U (nUT)l 2 (^,m,k)A k B m - (nUT)l 2 (k J <t,m)C k C m 


k=l 


m=l Y 


- Q 1 (mrr)l 4 (k,-t 5 m)B k C m + Q^kn) I 5 (k,£,m)B k Q 


In the derivation of Eqs . (19) through (21) a linear steady state velocity 
distribution, u = u e z, has been assumed, and the following definitions have 
been used: 


I 1 (k,-t') 


pl 

z sin(knz) cos (I-ttz) dz 

J 0 


( 22 ) 


I 2 (k,m,-t) 


pi 

sin(krrz) sin(nurz) cos(£ttz) dz 

J 0 


(23) 


I 3 (k,m,£) 


£ 


cos(krrz) cos(hutz) cos(£ttz) dz 


(24) 
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( 25 ) 


and 


I, (k,m,t) = J^z cos (krrz) sin(mnz) cos (-Dttz) dz 
4 0 


I (k,m,£) = J^z sin(knz) sin(mxrz) sin('trrz) dz 
5 0 


. v-1 _ 

Q-, = -TT~ U 
1 2y e 


q = - *-zi u 

^ 8y 2 e 


(26) 


The second order solutions are found by performing the operations 
defined by Eq. (18) . The following equations describe the behavior of the 
mode-amplitude function of the velocity potential: 


2 

d D. 9 cLD- r dD. dD„ 

= -(*>) ^ [_i . _i (t -T)] 


dt 


+ 2 j x feW 1 !! 1 ’ 4 ) - ¥ “e'- 1 ) 11 * 1 ] 3T 


N 


» r « an. 

+ E L(y-l)(mrr)^i (k,m , 1 ) — — D 


m=l 


dt m 


an -n 

2(kiT) (nm)l 2 (k,m,^) ^ dJ j 


(27) 


The space integrals defined in Eqs. (22) through (26) are evaluated 
numerically using a Simpson's rule algorithm 10 (see Appendix A) . The 
nonlinear behavior of axial mode instabilities are found by numeric all y 
integrating either Eqs. (19) through (21), or, in the case of moderate 
amplitude oscillations, Eqs. (27) . In order to carry out these computations, 
engine operating conditions (i.e., y, u , n and t) , and initial conditions 
must be specified. The behavior of the mode- amplitude functions is followed 
through the transient phase to the establishment of periodic oscillations. 

The perturbation flow field is then calculated using either Eqs. (10) through 
(12) or El* (13) . When Eqs. (27) are used to describe the unsteady flow, the 
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pressure perturbation at any location within the chamber is related to 9 hy 

5 6 

the following second order momentum equation : 


p'u,t) - 1 [? 2 (v0 - ^(v 2 *)] 


(28) 


A more detailed description of the solution technique outlined in the 
preceding discussion is included in Appendix B. lypical numerical solutions 
of these equations are presented and discussed in the following section. 

RESULTS AMD DISCUSSION 


Nonlinear Solutions 

Extensive computations have shown that the predicted nonlinear insta- 
bilities are dependent upon the engine operating conditions and independent 
of the nature of the initial disturbances. However, the computation time 
required to reach limit cycle conditions is reduced when the waveform of the 
initial disturbance is "close", in some sense, to the waveform of the limit 
cycle oscillation. For example, the computation time required to reach a 
discontinuous fundamental mode (1L) limit cycle oscillation is reduced when 
the assumed initial disturbance has a 1L discontinuous waveform. In this 
investigation it has been assumed that the engine is operating smoothly until 
t = 0, at which time a pressure disturbance is impulsively introduced inside 
the chamber. The velocity perturbation is assumed to be initially zero. Both 
spacially continuous and spacially discontinuous initial disturbances have 
been used. 

Typical transient and the resulting limit cycle oscillations are shown in 
Fig. 2. Here, initial continuous fundamental mode perturbations distort 
themselves into a discontinuous oscillation. When the amplitude of the initial 
disturbance is larger than the amplitude of the limit cycle oscillation, the 
transition to shock-wave-type behavior occurs within two cycles. On the other 
hand, when the initial amplitude is small a longer time period elapses before 
a shock wave is formed. In either case, the initial disturbances reach the 
same limit cycle conditions . These data were generated by solutions of the 
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Figure 2. Typical Transient and Limit Cycle Pressure Oscillations 
for a Given Engine Operating Conditions 


second order wave equation (i.e., Eq. (27)); however, the same behavior is 

3 

exhibited by the solutions of Eqs. (19) through (20) . It can also be shown 
that during the transient period the waves can change from one mode of 
oscillation to another. Consequently, if for given engine operating conditions 
the initial disturbance is not in the "proper" mode, then the solutions will 
adapt themselves to the operating conditions and the final periodic solution 
will be in the mode of oscillation that is unstable. In other words, no a 
priori knowledge of the behavior of the resulting instability is required in 
order to use the analytical technique developed in this report. 

The pressure envelope of the combustion instability oscillations is also 
defined in Fig. 2. The pressure envelope is simply the band of the peak-to- 
peak pressure amplitudes of the oscillations. The temporal behavior of the 
pressure envelope will subsequently be used to investigate engine triggering. 

It has been found in the course of this investigation that while the 

second order theory is capable of predicting the behavior of the final 

instabilities in linearly unstable engines , it is unable to predict the 

engine’s triggering limits. It is shown elsewhere^ that this difficulty is 

related to the mathematical structure of the resulting second order equations 

for the mode amplitudes. In view of these results, it was decided to use the 

second order theory, that requires considerably less computation time, to 

investigate the behavior of stable limit cycles in linearly unstable engines, 

while the third order theory will be used to study the behavior of triggering 

limits. A justification for this approach is presented in Fig. 3 where 

predictions of the second order and third order theories, for limit cycle 

peak-to-peak pressure amplitudes, are compared. It is shown in Fig. 3 that 

the predictions of both theories are in fair agreement over a wide region of 

peak-to-peak pressure amplitudes. A possible reason for the observed 

discrepancies is the different treatment of the nozzle boundary condition in 
o 3 

the two theories . It has also been shown^ that the waveforms predicted by 
the two theories are in good agreement. 

To determine the engine triggering limits, the minimum value of an initial 
disturbance required to initiate instability in a linearly stable region was 
determined numerically. For operating conditions where no disturbance can 
cause instability, the engine is said to be absolutely stable. Due to the 
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Second Order Solutions 



Figure 3. A Comparison of Second and Third Order Solutions 



above-mentioned shortcomings of the second order theory, the third order 
analysis is used herein to investigate the behavior of the triggering limits; 
the results of this investigation will determine the manner in which the 
problem’s nonlinearities modify the engine linear stability limits. The 
behavior of oscillations near a triggering limit are shown in Fig. 4. The 
upper plot in Fig. 4 was obtained for an initial disturbance larger than the 
triggering limit; it shows the growth of a pressure envelope. The lower plot 
in Fig. 4 shows the behavior of an initial disturbance whose magnitude is 
smaller than the triggering limit; the plot shows the decay of a pressure 
envelope. The plots presented in Fig. 4 indicate that the threshold 
disturbance amplitude required to trigger a linearly stable motor, at the 
operating conditions in question, lies between the amplitudes of these two 
initial disturbances . The threshold amplitude can be found by requiring a 
zero growth rate of the threshold disturbance and linearly interpolating the 
data shown in Fig. 4. 

The nonlinear behavior of fundamental mode instabilities can be summarized 
in an amplitude map of the type shown in Fig. 5 • This figure shows linear and 
nonlinear stability limits, and lines of constant peak-to-peak pressure 
amplitude. According to Fig. 5, triggering can be obtained in the narrow 
region between the linear stability limit (solid line) and the nonlinear 
stability limit (broken line) . The small extent of the triggering region is 
evident at f = 1.623 . Here, the vertical displacement of the nonlinear limit 
from the linear stability limit is only 6^ = n - n^, = - .02. The triggering 
region for above resonant conditions (i.e., t < 1) is also narrow, and is 
terminated at f = 2/3 where the second longitudinal mode becomes linearly 
unstable. It can be shown 3 that the concept of triggering becomes meaningless 
in a region where one of the modes present in the series expansion is 
linearly unstable. The significance of the parameter t Q /T, also shown in 
# j will be discussed m the following section. 

The objective of the preceding discussions was to provide an indication 
of the type of data that can be generated by the solution technique developed 
in this report. Detailed presentations of these and related studies can be 
found in Ref. 3* 
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Correlation with Experimental Data 


It has been shown in Reference 3 that the shape of the pressure waveforms 

depends upon engine operating conditions. Specifically, the pressure waveforms 

are dependent upon the proximity of the operating point to engine resonant 

conditions (i.e., to n . and t . ) . The observed behavior of the stable 

min min 

limit cycle pressure oscillations can be used to correlate the analytical 
results with experimental data. To accomplish this task, two waveform 
parameters are defined in Fig. 6. In this figure, the solid line shows the 
numerically computed pressure waveform, and the broken line is the mean pres- 
sure waveform used to determine the correlation parameters Ap ' (z ) and 

max r 

t /T(zJ; z r is the normalized axial location for which experimental pressure 
data is available. 

Once z r is specified, the analytical solution technique can be used to 
determine both the limit-cycle amplitude map and the dependence of t /T on t. 
Typical results are shown in Fig. 5- The values of Ap ' (z ) and t /t(z ) 
found from experimental pressure data are then used, in conjunction with the 
data presented in Fig. 5, to determine the engine operating conditions in 
terms of n and f. 

Semi-Empirical Pressure Waveforms 

A semi-empirical method for predicting the pressure waveforms has been 
developed. The objective of the semi -empirical method is to provide design 
engineers with a straightforward technique, requiring relatively little 
computation time, for predicting the nonlinear pressure waveforms. The semi- 
empirical correlation method is based on the observation 11 that the velocity 
potential, 9? can be approximated, at least for resonant oscillations, by the 
following series expansion: 


N 

9 = A, H k cos(kau,t) cos(krrz) 
k=l 


-O' 


( 29 ) 


where A^, a, and are found from computer-generated data. The nonlinear 
pressure waveform is then found from Eq. (28) . 
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The parameters A^, o/, and u}^ are found from the behavior of the mode- 
amplitude functions computed using Eq. (27) with a five term series expansion 
(i.e., N = 5 in Eq. (13)) • The parameters and cu are the maximum amplitude 
and the frequency of the fundamental harmonic, respectively. The exponent, a, 
accounts for the decrease in maximum amplitudes of the higher harmonics; it is 
found from an empirical log-log plot of maximum mode amplitude versus the mode 
number. 

Normally, a ten term series (i.e., N = 10, in Eq. ( 13 ) ) is required to 

q 

adequately predict the discontinuous waveforms . The required computation 
time is approximately proportional to the square of the number of terms 
retained in the series expansion. Consequently, the computation time required 
for the semi -empirical method is considerably .shorter than that required to 
solve directly for the pressure waveforms using the series solutions containing 
the unknown time -dependent mode amplitudes. 

Semi -empirical pressure waveforms are compared with computer generated 
solutions in Fig. 7- Ten terms were retained in Eq. (28) in the computation 
of the semi -empirical waveforms (i.e., N = 10 in Eq. (28)) . It is evident 
from the data shown in Fig. 7 that the semi -empirical method fails to 
reproduce the waveforms at off resonant oscillations. The probable reasons 
for this failure are: 

1. There is a slight phase shift between the various modes at off-resonant 
conditions . 

2. For off-resonant oscillations, the higher harmonics are both frequency and 
amplitude modulated. 

3 . For off-resonant oscillations, the higher harmonics may not obey the 
amplitude power law found by considering the behavior of the first few 
mode- amplitude functions . 


DISCUSSION AND CONCLUDING REMARKS 

An analytical technique has been developed for the analysis of nonlinear 
longitudinal combustion instabilities in liquid propellant rocket motors. The 
technique requires relatively little computation time and provides considerable 
insight into the physics of the problem. The method does not require any 
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Figure 7. Comparison of Numerical and Semi -Empirical Pressure Waveforms 



a priori knowledge of the final form of the instability. The method can 
predict triggering limits as well as the transient and final periodic 
behavior of the instability. Results predicted with the aid of this method 
agree with available experimental data. 

Results obtained with a second and a third order analyses show that the 
second order analyses describes the behavior of longitudinal combustion 
instability over a broad range of engine operating conditions. The third 
order theory showed that longitudinal instabilities can be triggered in a 
very narrow region outside the linear stability limits. The extreme 
narrowness of the nonlinearly unstable region suggests that from a practical 
point of view, the longitudinal stability limits of most engines are adequately 
described by the linear stability limits. 

A correlation technique, that can be used to correlate the analytical 
results with experimental data, and a semi-empirical method for predicting 
the waveforms of the instability, have been developed. 


26 



APPENDIX A 


PROGRAM SPAINT: EVALUATES THE SPACE INTEGRALS RESULTING 

FROM THE APPLICATION OF THE GALERKIN METHOD 

Statement of the Problem 

Program SPAINT uses a Simpson's rule integration algorithm to 
evaluate the space integrals resulting from the application of the 
Galerkin method. A linear ramp Mach number distribution, u(z) = u e X z, 
is used. The computed integrals are stored in a data file which is 
used as input data in Program WAVES. The program user must specify 
the step size to be used in the numerical integration, and the number 
of terms retained in the series expansion(s) of the dependent variable(s) . 

The space integrals to be evaluated are given in Eqs. (22) 
through (26) . The following definitions are made for the purpose 
of computer storage assignment: 


Array 

Integral 

pi 

Index (K) 


T2(1,N,L) 

= z sin(nrrz) cos(Lttz) dz 

tJ 

0 

pi 

0 

(A.l) 

T3(1,N,M,L) 

= sin(nTTz) sin(mrrz) cos (£ttz) dz 

u 

0 

pi 

1 

(A. 2) 

T3(2,N,M,L) 

= cos (mrz) cos (mrrz) cos (Lttz) dz 

2 

(a.3) 


0 
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Array 


Integral 


Index (K) 


T3(3,N,M,L) = Jz cos(nnz) sin(mrrz) cos(4ttz) dz 3 (A. 4) 


0 


r 1 

T3(4,N,M,L) = z sin(nnz) sin(mrrz) sin('t/Trz) dz 4 

0 


(A.5) 


The array indices N, M, and L vary from one to NEQ, where NEQ 
is the number of terms retained in the series expansion(s) of the 
dependent variable(s) . As coded in this report, NEQ ^ 10. It is 
recommended that a value of NEO, = 10 be used regardless of the n umb er 
of terms in the series. Ihe reason for this choice is discussed in 
the section of this appendix entitled "Recommendations on Program 
Usage". 

A standard Simpson’s rule numerical integration algorithm (see, 
for example, Conte 10 ) is used to evaluate the integrals. In this 
procedure, the interval [o,l] is divided into 2N subintervals of 
length h and the integral is evaluated using the following equation: 


. f W 4* - | [f 0 + 4f 1 + 2f 2 + 4f 3 + . . . + 4^ + f 2N ] 

0 


The error involved in this numerical integration scheme is of the 

4 

order of h . Ihe user specifies h, and h must be such that the interval 
[0,1] is divided into an even number of subintervals. 
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Input and Output 

The required input data consist of the number of terms in the 
series expansion(s) of the dependent variable ( s) , NEQ, and the integra- 
tion step size, HI. The input data is read into the computer from 
two data cards : 

Card 1: NEQ, integer, is right justified in columns 1-10 (Format 110) 

and NEQ ^ 10 

Card 2: HI, floating point number, in columns 1-10 (Format F10.0) 

The computed integrals are stored in an assigned data file (see 
the section on the Deck set-up) and are printed in a straightforward 
output format. The notation used in the printed output is self- 
explanatory: L, N, and M are array indices (M = 0 for integral (A.l)) 

and K is the index which defines the integrand (e.g., K = 0 for 
integral (A.l) , etc.) . 

A typical set of input data and a portion of the printed 
output axe respectively shown in Tables A.l and A. 2. 

Deck Set-up 

The deck set-up described herein is for the Univac 1108 Exec 8 
system used at Georgia Tech. The manner in which data files are 
assigned might be different at other computer facilities. The important 
thing to note is that the data file number (i/O unit) assigned to the 
output data of this program is used as the input data file number in 
program WAVES. This program uses i/O unit 2 to store the data file. 

Deck Set-up: 

1. Run Card (i.D. Card) 
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TABLE A. 2. Sample Output From Program SPAINT 
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X 
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2. i/O unit assignment cards. 


3- Main Program, MAIN. This program reads the input, calls subroutine 
SUMM, and outputs the computed integrals . 

4. Subroutine SIJMM. This program specifies the integrand function, 
f (x) , and calls subroutine SIMPSN. 

5* Subroutine SIMPSN. This program performs the Simpson rule integra- 
tion of f(x) . f (x) is defined in the External Real Function 
Subprogram FOFX. 

6. Real Function Subprogram FOFX. This program defines the integrand 
function f(x) according to the integral index, K. 

7 . Input Data Cards . 

Recommendations on Program Usage 
Experience with this program has shown that an integration step 
size of HI = .02 produces good results. Although NEQ can be varied 
from 1 to 10, it is recommended that NEQ = 10 be used for the following 
reason: Using this approach, one data set can be used to compute 

nonlinear solutions (using program WAVES) for values of NEQ between 
one and ten. Program WAVES is set-up to use the output generated by- 
program SPAINT in this manner. In summary, it is recommended that 
values of HI = .02 and NEQ = 10 be used. Approximately 60 seconds of 
computation time on a U-1108 are required in this case. 
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oooonooooooooooooooooo 


FORTRAN Listing of Program SPAINT 


THIS PROGRAM EVALUATES THE INTEGRAL OF F(X) FROM 0 TO 1 
USING SIMPSON RULE 

THE MAIN PROGRAM READS THE INPUT* CALLS SUBROUTINE SUMM* 

AND OUTPUTS THE COMPUTED INTEGRALS. THE INTEGRALS ARE 
PRINTED AND STORED IN FILE 3 USIN6 THE FASTRAN SYSTEM. 

THE F ( X ) ARE DEFINED WITH THE PRINTED OUTPUT 
THE SIMPSON RULE INTEGRATION IS PERFORMED IN 
SUBROUTINE SUMM. THE F < X ) ARE DEFINED IN THE EXTERNAL 
FUNCTION SUBPROGRAM FOFX. 

INPUT DATA 

CARD 1 IN COL. 1-10 THE STEP SIZE* HI (ABOUT .01 TO .02) 
CARD 2 RIGHT JUSTIFIED IN COL. 1-10 THE NUMBER OF 
TERMS IN THE SERIES EXPANSION NEQ< OR = 10 

THE OUTPUT DATA IS DEFINED IN THE PRINTED OUTPUT 

THE COMPUTATION TIME ON THE U-1108 IS ABOUT 60 SEC FOR 
HI = .02 AND NEQ = 10. 

COMMON/ INTER/ T2 ( 1 * 10 . 10 ) * T3 < 4 * 10 . 10 * 10 ) 

400 FORMAT (0110) 

402 FORMAT (10X* 'OUTPUT FORMAT INTEGRAL FROM 0 TO 1 OF F(X)'*/ ) 

403 FORMAT <10X*'K=0 IS F<X) = X*SIN(N*PI*X)*CoS(L+PI*X) ' ) 

405 FORMAT (10X*'K=1 IS F(X) = SIN(N*PI*X)*SIN(M*PI*X)*COS(L*PI*X) ' ) 

406 FORMAT (10X.'K=2 IS F<X) = COS(N*PI*X> *COS(M*PI*X> *COS(.L*PI*X> • > 

407 FORMAT (10X*'K=3 IS F(X) = X*COS ( N*PI *X ) *SIN ( M*PI *X > *COS ( L*PI*X ) • ) 

408 FORMAT <10X*'K=4 IS F(X) = X*SIN ( N*PI *X ) *SIN ( M*PI *X > *SIN ( L*PI *X > • ) 

410 FORMAT (8FI0.0) 

430 FORMAT (1H1. lOX.'SPACE INTEGRALS STEP SIZE = '*F5.3* 

1 4X* 'L =' * 12*/) 

440 FORMAT (215* 10E10.4) 

450 FORMAT (/*0X*2HN=* 15* 9110* ) 

460 FORMAT (' K M'/) 

800 FORMAT (5E15.8) 

READ (5*410) HI 
READ (5*400) NEQ 
C INTEGRATION OF SPACE INTEGRALS 

CALL SUMM (NEQ* HI) 

DO 200 L=1 * NEQ 
WRITE (6*430) HI * L 
WRITE (6*402) 

WRITE (6*403) 

WRITE (6*405) 

WRITE (6*406) 

WRITE (6*407) 

WRITE (6.408) 

WRITE (6*450) ( I * 1 = 1 * NEQ) 

WRITE (6*460) 

M = 0 
K = 1 
J = 0 

WRITE (2.800) (T2(K*N*L) *N=1*NEQ) 
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WRITE 

( 6 1 

440) 

JrM» 

<T2< 

K#N#L) 

» N= 1 * 

NEQ) 


DO 220 

K: 

= 1*4 







DO 230 

M: 

:l»NEQ 







WRITE 

( 2 1 

>800) 


( T3 ( 

K # N # M » 

L> *N= 

1»NEQ) 

230 

WRITE 

( 6 1 

>440) 

K. M» 

<T3< 

K » N » M » 

L> »N= 

1 » NEQ ) 

220 

CONTINUE 







200 

CONTINUE 








WRITE 

( 2 1 

>800) 

HI 






STOP 









END 









SUBROUTINE SUMM (NEQt HI) 

COMMON/ INTER/ T2 ( 1 , 10 . 10 ) » T3< 4 > 1 0 # 10 • 10 ) 
NSM=1./HI + 1.01 
00 100 K=l#5 
IF (K.GT.l) GO TO 1 
MSTP=1 
GO TO 3 
1 MSTP=NE0 

3 00 200 L=1 • NEO 
AL = L*3. 14159 
DO 210 N=1#NEQ 
AN = N*3. 14159 
DO 220 M=1 • MSTP 
AM = M*3. 14159 

CALL SIMPSN (K»AL*AM» AN*5UM* NSMrHI) 

IF (K.GT.l) GO TO 4 
T2(KrN»L) = SUM 
GO TO 5 

4 KK =' K-l 

T3 ( KK » N * M* L ) = SUM 

5 CONTINUE 
220 CONTINUE 
210 CONTINUE 
200 CONTINUE 
100 CONTINUE 

RETURN 

END 



- NSM*HI> 


SUBROUTINE SIMPSN ( K r AL » Ay # AN » SUM t 
EXTERNAL FOFX 
X = 0.0 
SUM =0.0 
DO 1 I = 1 » NSM 
C = 1.0 

IF (I.EO.l) GO TO 2 
IF (I.EQ.NSM) GO TO 2 
C = 4.0 

ID - 2* ( 1/2 ) - I 
IF (IO.EO.O) GO TO 2 
C = 2.0 

2 SUM = SUM + C*FOFX<K»X»AL»AM.AN> 

1 X = X+Hl 

SUM = HI *SUM/3 • 0 

RETURN 

END 


1 

2 

3 

4 

5 

100 


REAL FUNCTION FOFX ( K t X » AL * AM> AN ) 
GO TO ( 1 r 2 * 3 f 4 * 5 ) » K 


FOFX = X*SIN(AN+ 

GO TO 100 

FOFX = SI N ( AN*X ) 

GO TO 100 

FOFX = COS ( AN*X ) 

GO TO 100 

FOFX = COS ( AN*X ) 

GO TO 100 

FOFX = SI N ( AN*X ) 

CONTINUE 

RETURN 

END 


X ) *COS ( AL*X ) 

*SlN( AM*X)*COS( AL*X> 
*C05( AM*X) *C0S( AL*X> 
*SIn(Avi*X)*COS( al*x)*x 
*SIN< AM*X) *5IN( AL*X)*X 
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APPENDIX B 


PROGRAM WAVES: COMPUTES THE COMBUSTION 

INSTABILITY OSCILLATION WAVEFORMS 

Statement of the Problem 

Program WAVES computes the combustion instability oscillation 
■waveforms for combustors having a linear steady state velocity distri- 
bution, u(z) = u g z, for which u e is small. Before this program can be 
used, the space integrals must be evaluated using program SPAINT. The 
computed integrals , together with the specification of the engine 
operating conditions (i.e., n, t, u g , y, etc.), initial conditions, 
and certain program control numbers, make up the required input data 
for program WAVES. 

Program WAVES performs the following functions : 

1. For an initial peak pressure amplitude, initial values of the mode- 
amplitude functions are computed. 

2. The time-dependent mode- amplitude functions are found by a Runge- 
Kutta-type numerical integration. 

3. Perturbation pressures and velocities are computed. 

4. A check for limit cycle conditions is made. 

5. Printed and/or plotted output data is generated. 

The program provides the user with various options. For instance, 
function ( 3 ) may be omitted if only the behavior of the mode- amplitude 
functions is desired. Similarly, function (4) is omitted when only the 
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transient behavior of the instabilities is required. The use of 
these and other user options are discussed in this appendix. 

Two nonlinear solutions have been developed in this report: 

(l) a second order analysis using a nonlinear wave equation, and (2) a 
large amplitude analysis using a set of three conservation equations. 
Consequently, two computer programs are required. These programs 
have been written in a manner which permits a good deal of commonality. 
In particular, the required input data is the same for all programs. 

In order to achieve the commonality between the programs , the 
definitions shown in Table B.l have been made. 

The relations defining the behavior of the functions A (t) , 

B (t) , and C (t) are listed in Table B.2. 

Program WAVES consist of 11 elements: MAIN, START, POFX, TREND, 

FLOW, P0UT2, POUT, RUNG, EQTN, PRMTRS , and W0UT1. The first seven 
elements are the same for the three nonlinear solutions. The last four 
elements are different for each nonlinear solution technique. The 
functions performed by these elements are discussed in the following 
paragraphs . 

MAIN: Element MAIN serves the twofold functions of (l) reading 

the data required to compute the nonlinear waveforms, and (2) calling 
the required subroutines . 

START and POFX : These two subroutines provide the initial values 

of the mode-amplitude functions required for the integration of the 
ordinary differential equations describing the behavior of the mode- 
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TABLE B.l. Definition of the Mode-Amplitude Functions 
Used in Program WAVES 


Array Parameter 


Description 


a(n) 

B(N) 


C(N) 


A n (t) Specific volume mode-amplitude function, used 

only in the large amplitude analysis . 

B (t) Pressure mode-amplitude function. In the analysis 

using the nonlinear wave equation, B (t) represents 
the time derivative of the velocity potential- 
mode- amplitude function. 

C (t) Acoustic-type velocity mode -amplitude function. 

In the nonlinear wave equation solutions, C (t) 
represents the velocity potential mode- amplitude 
function. 


TABLE B.2. Equations Governing the Mode-Amplitude Functions 


Parameter Equation Number 

Wave Equation Third Order Equations 


A (t) 

— 

19 

B (t) 
n v ' 

27 

20 

C (t) 
n ' 

27 

21 
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amplitude functions. It is assumed that the conibustor is operating 
in a steady manner until time t = 0, at which time a pressure disturb- 
ance is impulsively introduced inside the combustor. The perturbation 
velocity at t = 0 is zero. The user may specify a spacially continuous 
initial pressure disturbance in any axial mode, or a spacially discon- 
tinuous fundamental mode disturbance, with the discontinuity located 
at z = .5 at t = 0. The analytical expressions used to find the initial- 
conditions, found by a Fourier analysis of the initial waveform, are 
given in the following equations : 

(1) Spacially Continuous Pulse in the -tth Axial Mode. 


c.(t = 0) = 0 


l = 1,...,N 


(B.l) 




(B.2) 


(2) Spacially Discontinuous Pulse. 


C^(t = 0) = 0 


l = 1, . . . ,N 


(B.3) 


B t (t * °> = psj sin (r) 


(B.4) 


where in both cases, 


C p ( t) = B » ( t) — 0 , for -T ^ t <0, •£' — 1,...,N 


(B.5) 
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An initial condition on A^(t) is required for the large amplitude 
analysis. Consideration of the linear behavior of the system shows 
that v = -p/y; consequently the following relation is used for an 
initial condition on A^(t = 0) : 

A t (t = 0) = -B^t = 0)/y (B.6) 

In the solution using the nonlinear wave equation, Eqs. (B.2) and 

(B.4) merely approximate the spacial dependence of the initial impulse. 

Specifically, these equations are based on a linear representation of 

the initial disturbance, and the computed wave amplitude differs by a 

factor of y from the specified p^. 

TREND : This subroutine determines whether or not limit cycle 

conditions have been reached. This task is accomplished by evaluating 

NEQ 

the summation S = B n (t) and examining the behavior of the summation. 
Note that S represents the behavior of the injector face pressure. 
Subroutine TREND performs the following functions : 

1) Determines the maximum (positive) peak amplitude of one cycle of S. 

2) Finds two successive average values of S for two cycles, S-^ and S^, 
respectively. 

3) Compares the absolute difference, | AS | , between the two successive 
averages with a user specified percentage, e, of the latter value of the 
average S. If the |as| < eS^ then limit cycle conditions have been 
reached . 

4) Makes the appropriate change in the internal program control index 
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which tells the program that limit cycle conditions have been 
reached. 

FLOW: Subroutine FLOW computes the summations used to find the 

perturbation flow field, outputs the computed pressure and velocity, 
and calls subroutine P0UT2. The summations computed are: 

NEQ 

SUMA = E A (t)cos(nnz) 
n=l n 


SUMS 


NEQ, 

S 

n=l 


B (t) cos (nrrz) 


NEQ 

SUMO = E 
n=l 


C n (t) sin(nrrz) 


SUMO = L (nn)c (t) sin(nrrz) 
n=l n 

Fnese summations axe used in subroutine PRMTRS to calculate the 
perturbation flow field. 

P0UT2 and POUT : Subroutines POUT and P0UT2 plot the temporal 

behavior of B(N) (the pressure mode-amplitude functions) and the 
temporal behavior of the pressure oscillations, respectively. The mode- 
amplitude functions to be plotted are specified by the user. The axial 
location (s) of the pressure plots are also user specified. The programs 
have been developed for use on a CALCOMP plotter. 
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RUNG: Subroutine RUNG is a modified Runge-Kutta numerical 


integration algorithm. The modification accounts for the presence 
of the retarded time variable. In this modification, the retarded 
variable is treated as a known quantity; that is, it is treated in 
the same manner as the independent variable. Two algorithms, based 
on the Runge-Kutta equations developed in reference (12) , are used. 
One algorithm is used to integrate a set of second order O.D.E.'s; 
the other is used to integrate a set of first order O.D.E.'s. The 
required expressions are given in the following equations : 

(l) First order O.D.E.'s; 

y[3 +1) = y<"» + | {K^ ♦ ^ J 


where 



\ / 2 )- y n (t-T*/ 2 )] 

t n 

k = / 2 )> y n (t-^/ 2 )] 

t n 

\ = “t’Kvs )’ y^- 7 *)] 

Aj n 
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and where 


yf = y t (t) 


= y^ft-th) 


and is the function evaluated at t. 

(2) Second order O.D.E.'s; y^' = y n^ t- ^] 


y t 


i (J+1) - yl (d) + 5 {\ + \ + i\ * \]} 


y(J +1) = y<J> ♦ h{y^> + | [k + K + K ]} 


l t D l~ 


where 



y n’ y n 





, h / 
+ 2 y n 




T + 




, h / 
+ 2 y n 




T + 


MB' MB' 





+ hy. 


n 



+ ). - t + h )] 
J n 


axid where 


y[ 3) = y t (t) ; y^> = y'(t) 

yp +1 ^ = y t (t-Hi) ; y^ +1 ^ = y|(t-Hi) 

The equations defining the numerical integration of a set of 
first order O.D.E.'s is used in the solutions of the conservation 
equations. The second order O.D.E. equations are used to solve the 
nonlinear wave equation. The functional form of f^ is defined in 
element EQTN. 

In order to use the equations with the retarded variable, the 
integration step size, h, must be selected such that h divides the 
time lag, f, into K equal increments. Thus f = Kh, and the retarded 
variables become: 


y n (t - n 


y (t 

J n 


Kh) 


y n (t - r + §) = y n (t 


Kh +| ) 


y n (t - T + h) = y n (t - Kh + h) 
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It has been found that an integration step size of the order 
h — .05 produces satisfactory results. The program selects the 
integration step size hy forming the ratio t/. 05 > rounding off the 
result to the nearest integer, and dividing f hy the resulting integer, 
that is : 


integer = (f/. 05 ) + *01 
h = f/integer. 

The computation of h is performed in element MAIN. 

EQTN: Subroutine EQTN defines the functions, f^, used in 

subroutine RUNG to evaluate the K terms . The particular equations 
defined in EQTN depend upon the problem under consideration (i.e., 
nonlinear wave equation, etc.) . These functions are defined in Table 
B.2. 

PRMTRS : Subroutine PRMTRS uses the summations, SUMA, SUMB, 

SUMC, and SUMU, computed in subroutine FLOW to calculate the perturbed 
flow field. The current program is coded to compute the perturbation 
pressure and velocity, using the following equations: 

(l) Nonlinear wave equation solutions: 

u'(z,t) = -SUMU 

p'(z,t) =| [sUMB(SUMB-2) + SUMU(2u(z) - SUMU) ] 
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(2) Second order conservation equation solutions: 

p 7 ( z j t) = SUMB 


u'(z,t) = SUMC + ^ u Z SUMB 
v 5 ' 2y e 


(3) Large amplitude solutions: 

p'(z,t) = SUMB 


u / ( z j t) = SUMC + [l - SUMB] ^ U e Z SUMB 

WOUT1: This program writes the output of the mode -amplitude 

functions . 


Input Data 

The required input data consist of the integral values computed 
by program SPAINT, the engine operating conditions, and program control 
numbers. The data from program S PAINT is automatically read from data 
file 2. The remaining data is read from user supplied data cards. 

These cards are described in this section. 
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Card 1 (Format 8H0) 


Column" 1 " 

Term 

2 

Data Type 

Information 

Restrictions 

10 

NEQ 

I 

No. of terms in the 
series expansion of the 
dependent variables 

£10 

20 

NX 

I 

No. of axial locations 
at which flow field is 
to he computed 

£11 

30 

LIN 

I 

LIN = 1 to compute 
linear solutions 
LIN ^ 1 nonlinear 
solutions 


4o 

IP LOT 

I 

IPL0T = 1 if any data 
is to be plotted 
IPL0T / 1 no plots 


50 

Card 2 

INPT I 

(Format 8110) 

INPT = 1 to write the 
space integrals 
INPT £ 1 space integrals 
are not written 


Column 

Term 

Data Type 

Information 

Restrictions 

10 

LC1 

I 

LC1 = 1 to write the 
mode- amplitude functions 
LCl ^ 1 mode -amplitude 
functions are not written 


20 

LC2 

I 

LC2 = 1 to plot pressure 
mode- amplitude functions 
LC2 = 4 no plot of mode- 
amplitudes 


30 

lc4 

I 

Number of terms to be 
plotted 

£10 

40 

LC5 

I 

Incremental index between 
terms to be plotted 

^ 9 

1. For 

integer 

data, indicates 

the column in which data is 

right 


justified. 

2. I denotes integer data; F denotes floating point (decimal) data. 
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Card 3 (Format 8110) 


Column 

Term 

Data Type 

Information 

Restrictions 

10 

LP1 

I 

LP1 = 1 to calculate 
p and u 

LP1 ^ 1 flow field is 
not calculated 


20 

LP2 

I 

LP2 = 1 to write p ' 
and u 

LP2 / lp and u are 
not written 


30 

IP3 

I 

LP3 = 1 to plot p ' vs t 
LP3 = 4 no flow field 
(p ) plot 


40 

lp4 

I 

Number of axial locations 
at which p vs t is to 
plotted 

<;4 

Card 4 

(Format 8110) 



Column 

Term 

Data Type 

Information 

Restrictions 

10 

NTAU 

I 

Number of f to he run 


Card 5 

(Format 8F10 .0) 



Column 

Term 

Data Type 

Information 

Restrictions 

1-10 

UE 

F 

Exit Mach number 

small, « 1 

11-20 

GAMMA 

F 

Specific heat ratio 


21-30 

EPS 

F 

Limit cycle amplitude 
percent error 

EPS = er( .01) 

Card 6 

(Format 8F10 .0) 



Column 

Term 

Data Type 

Information 

Restrictions 

1-10 

TBEGIN 

F 

Normalized time at which 
output is begun, and at 

see 

discussion 


which flow field calcula- 
tion is started 
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Column 

Term 

Data Type 

Information 

Restrictions 

11-20 

TEND 

F 

Normalized time at which 
computations are 
terminated 

see 

discussion 

21-30 

TIMCY 

F 

Normalized time at which 
limit cycle check is 
begun 

see 

discussion 

31-40 

DELTAT 

F 

Normalized time increment 
for output of limit cycle 
conditions 

see 

discussion 

41-50 

TSMP 

F 

Normalized time at which 
plot of pressure mode- 
amplitude is begun 


51-60 

DELPT 

F 

Normalized time increment 
for plot of pressure mode- 
amplitude, B(N) vs t 

see 

- discussion 


Discussion of Card 6: 

(1) BEGIN must be greater than or equal zero. TEND must be such 
that the ratio (TEND-TBEGIN) /H is less than 300. This ratio can he 
estimated using a value of H = .05- Experience has shown that a time 
increment of TEND-TBEGIN * 12. is sufficient to determine the behavior 
of the solutions . 

(2) If a limit cycle check is not desired, then set TLYMCY > 

TEND. 

( 3 ) DELTAT must he such that DELTAT/H < 300. Usually, a 
DELTAT = 6 is sufficient to verify that limit cycle conditions have 
been reached. In this case, approximately three fundamental mode cycles 

are computed. 

(4) If a limit cycle check is made, and if limit cycle condi- 
tions are found, TSMP is automatically set equal to the initial time at 
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which limit cycle conditions are found, if LC2 = 1. 

(5) DEEPT must he such that DEEPT/H < IDO. Good results have 
been obtained using DEEPT = 3 .9. 

(6) If a limit cycle check is made, and limit cycle conditions 


are not found, the data output begins at TBEGIN and ends at TEND. 
Card 7 (Format 8F10.0) 

Column 

Term Data Type 

Information 

Restrictions 

1-10 

X(l) F 

Axial location at which 
p and u are computed 

£1 

11-20 

X(2) F 

Axial location at which 

SI 

• 

• 

• • 

• • 

p and u 7 are computed 

• 

• 

• 

• • • 

: x(nx) : 

Discussion: 

If NX > 8, then two cards 

• 

• 

• 

• 

: are required to complete 

• 

• 

• 

• 

the input of 

X(I) . 
Card 8 

In this case, X(9) is in 

(Format 8110) 

This card is included in 

columns 1-10 of card 7B, and so on. 
the data set only when EP3 = 1. 

Column 

Term Data Type 

Information 

Restrictions 

10 

IPX(l) I 

Index of X(l) at which 
a p vs t plot is made 

£10 

20 

• 

IPX (2) I 

Index of X(l) at which 
a p vs t plot is made 

£10 

• 

ko 

IPX(EP4) I 

Index of X(l) at which 
a p vs t plot is made 

£10 
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Discussion: 


Plots can be made at any four (or fewer) axial locations at which 
p ' is calculated. 

Card 9 ( 8) 1 (FID.O, 2110) 


Column 

Term 

Data Type 

Information 

Restrictions 

1-10 

TAU 

F 

Sensitive time lag 


20 

NNB 

I 

Number of n to be run 
at the specified t 

<■10 

30 

LCUT 

I 

Highest mode in which 
energy feedback is 
permitted 

see 

discussion 


Discussion: 

This number is used to eliminate the secondary zones of insta- 
bility. For fundamental mode investigations, LCUT = 2 is usually 
appropriate. For t > 1, energy feedback is only permitted in the 
fundamental mode. 


Card 10 

(9) (8F10.0) 


Column 

Term 

Data Type 

Information 

1-10 

ANR(l) 

F 

First value of n 

11-20 

ANR(2) 

F 

Second n 


ANR(NNB) F Final value of n 


Restrictions 


Discussion: 

If NNB > 8 , then two cards are used to input the ANR(l) . 

Number in parenthesis is the card number if card 8 (iPX(l) card) is 
omitted. 
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Card 11 

( 10 ) ( 8110 ) 




Column 

Term 

Data Type 

Information 

Restrictions 

10 

NPI 

I 

Number of initial 
disturbances for each 
n-f condition 


Card 12 

(11) (F10 .0 

, 8110 ) 



Column 

Term 

Data Type 

Information 

Restrictions 

1-10 

PI 

F 

Initial disturbance 
peak amplitude 


20 

IPOP 

I 

If IPOP ^ ID, then an 
initial disturbance in 



the IPOP mode is gener- 
ated 

If IPOP = 11, then a 
spacially discontinuous 
fundamental mode wave, 
with the discontinuity 
at z = . 5 , is generated 


This completes the description of the input data cards. If 
KPI > 1, then card 12(11) is repeated NPI times. When MB > 1, then 
cards 11(10) and 12(11) must he repeated MB times. Similarly, when 
NTAU > 1, card 10(9) through 12(11) must he repeated NTAU times. An 
example input data set is shown in Table B. 3 . 

Using the input data shown in Table B. 3 , program WAVES performs 
the following functions : 

1. Nonlinear solutions are found at two axial locations using eight 

term expansion(s) . The exit Mach number is u = 0.2, and v = 1.2. 

e 

2. The mode-amplitude functions are printed, and the first pressure 
mode-amplitude function is plotted. 

3. The perturbation pressure and velocity are computed at z = 0.0 and 
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COLUMN 


1-10 11-20 21-30 


0.2 

30.0 

0.0 

1.30 

1.18 

.025 

.05 

. 15 
1.0 
1.10 

. 1 


8 

1 

1 

2 


42.0 

0.25 


1 


2 


1.30 


1 


1 


2 

1 

1 


0.01 

5.0 


2 


1 

1 

11 

1 


11 




41-50 


51-60 


6.0 


30.0 


3.9 


2 



z =0.25. The results are printed and the temporal dependence of 
the pressure oscillation is plotted at z =0.0. 

4. A limit cycle check is initiated at t = 5*0. If limit cycle 
conditions are reached prior to t = 30 .0 , the required data is 
output in a time interval of At = 6.0 after the establishment of 
limit cycle conditions. On the other hand, if a limit cycle is not 
reached by t = 30.0 , the data is output in the time interval 

30 k t ^42.0. The pressure mode- amplitude function is plotted 
over at time interval of At = 3.9. 

5. Solutions are to be calculated for two values of t. At the first 
T (t = I.30) , the computations are to be made for two n (n = 1.18 
and n = I.30) . The computations at t = I.30, n = 1.18 are to be 
made using two initial disturbances; a .025 and a .05 peak 
amplitude 1L pressure wave. The computations at t = I.30, n = I.30 
are made for a discontinuous 1L pressure wave of peak amplitude 
equal to .15. At the second t (f = 1.0) , the computations are made 
for an n = 1.10 and a discontinuous, .1 peak amplitude pressure 
wave. 

6. In both cases, energy feedback is only permitted in the first two 
axial modes . 


Output Data 

The following data output options are available: 

(1) IEPT = 1 causes the space integrals used in the computations to be 
written. 

(2) LC1 = 1 results in a tabulated output of the mode-an^litude 
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functions . 


(3) LP2 = 1 results in the listing of p 7 and u* as functions of t at 
each axial location specified by X(l) . 

(4) LC2 = 1 causes plots of B(N) vs time to be made, with the N's 
specified by the user. 

(5) LP3 = 1 causes plots of p ' vs time to be made at the axial 
locations X(l) specified by the indices IPX(l) . 

The output limitations have been discussed in the data input section of 
this appendix. The output symbols are described in Table B.4. Portions 
of an example output is shown in Table B. 5 . 

Deck Set-up 

The data set described herein is for the Univac 1108 Exec. 8 
system as used at Georgia Tech. Tne important points are: 

1. Unit 2 must be assigned to the data file SPAINT. 

2. Unit 3 must be assigned to the CALCOMP PLOT subroutines. 

It is convenient to group the program elements in the sequence 
in which they are discussed in the first section of this appendix 
(i.e., page 37 ) . The program is then adapted to the solution of a 
particular formulation of the problem (i.e., second order wave equation, 
etc.) by changing the last four subroutines. 
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TABLE B.4. Output Symbols 


Symbol 

A 

AP 

B 

C 

LINEAR 

L 

M 

N 

NEQ 

P 

PINITIAL 

TAU 

Z 


(1) potential mode-amplitude function, or 

(2) specific volume mode- amplitude function 

time derivative of the potential mode- amplitude 

pressure mode -amplitude function 

velocity mode -amplitude function 

LINEAR = 1, solutions are linear 
LINEAR ^ 1, solutions are nonlinear 

axial mode number 

axial mode number 

(1) axial mode number, or 

(2) interaction index 

number of terms used in the solutions 
normalized perturbation pressure 
peak amplitude of the initial disturbance 
sensitive interaction index, t 
axial station 
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TABLE B.5 (cont.) . Sample Output From Program WAVES 

Initial Pressure Pulse 


INITIAL 

PRESSURE DISTRI3UTI0N 


Z 

P 

.00000 

.09216 

.10000 

• 1 0 6 8 4 

.20000 

.09628 

.30000 

.09790 

.40000 

.11386 

.50000 

.oonoo 

.60000 

-.11386 

*7-0 000 

-.-0-9790 

.80000 

-.09623 

.90000 

-• 10684 

1.00000 

-.09216 

1.10000 

--.10684 



TABLE B . 5 (cont.) . Sample Output From Program WAVES: 

Part of the Mode-Amplitude Output 
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TABLE B.5 (concluded) . Sample Output From Program WAVES: 

Part of the Perturbation Flow Output 


FLOW PRArAMEjERS 

Z= .000 


LINEAR:: 

2 

NEGi= 8 N= 1.10000 

TAU= l.OOOOn 

Exit maCh= .200 

SAWA =1.200 PTMITIAL =.1000 


TIVE 

PRESSURE 

VELOCITY 


-.000 

-.1055+0a 

.0000 


.050 

-.1104+00 

.0000 


.100 

Iln5+0 0 

• noon 


.150 

-.1185+00 

.0000 


. 20Q 

-« I0s2+0u_ 

.0000 


.250 

-.9664-01 

.0000 


. 300 

-,98fj7-0l 

.0000 


.350 

-.1142+00 

.0000 


.'400 

-.1213+00 

. nono 


.450 

-. 8656-01 

.0000 


.500 

-.2861-02 

.0000 


. 550 

.9295-01 

. 0000 


.600 .1479+00 

.0000 


. 650 

.1497+00 

.0000 


.700 

. 1144+00 

.0000 


,750 

.9950-01 

,0000 


.800 

.1125+00 

.0000 


.850 

, 13Q8+00 

.0000 


.900 

. 13o9+0n 

.0000 


.950 

.1168+00 

.0000 


1.000 

. 10°3+0n 

. 0000 


1.050 

,1143+00 

.0000 


1 . 10(1 

.1176+On 

.0000 


1.150 

, 1094+00 

. O000 


1 • 2 0 0 

.9693-01 

. 0000 


1.250 

,9826-01 

.0000 


1.300 

. 1 0 0 7 + 0 0 

.0000 


1.350 

.9221-01 

.0000 


1.400 

.7471-01 .0000 


1.450 

.5890-01 

.0000 


1.500 

»_3&i_9-0l 

.0000 


1.550 

-.4975-02 

.0000 


1.500 

-.6246-01 

.0000 


1.650 

-.1002+00 

.0000 


1.700 

-.1010+00 

.0000 


1.750 

-.3766-01 

.0000 


1.800 

- , 9498-0 1 

,0000 

• 

1.850 

-,12ol+Oo 

.0000 


1.900 

-.1309+00 

. (toot) 


1 « 950 

-. 1133+00 

,0000 


2,000 

-.9116-01 

,0000 


2.050 

-.9565-01 

.0000 




oooooooooooooooooooooo 


FORTRAN Listing of Program WAVES 


THE 

1, 


THE SPACE INTEGRALS ARE STORED IN THE ARRAYS TJ> AND T3. PROVISION 
IS MADE FOR ONE N BY N INTEGRAL* AND FOUR N BY N BY N INTEGRAL • 

MORE INTEGRALS CAN BE TREATED BY CHANGING THE APPROPRIATE DIMENSION 

STATEMENT. 

THE MODE AMPLITUDES ARE STORED IN THE ARRAYS A* B* AND C. THE RETARDED 
VARIABLE IS STORED IN ARRAY BS. THE RETARDED VARIABLES REQUIRED AT 
THE 1 INTEGRATION STEP IN QUESTION ARE STORED IN ARRAYS BR1 * BR2* AND 3R3. 
TERMS STORED IN THESE ARRAYS DEPENDS ON THE PROBLEM FORMULATION. 

, FOR THE NONLINEAR WAVE EQUATION 
BLANK 

TIME DERIVATIVE OF MODE AMPLITUDE 
MODE-AMPLITUDE FUNCTION 

THE SECOND ORDER CONSERVATION EQUATIONS 
3LANK 

PRESSURE MODE-AMPLITUDE 
VELOCITY MODE-AMPLITUDE 
THE LARGE AMPLITUDE ANALYSIS 
SPECIFIC VOLUME MODE-AMPLITUDE 
PRESSURE MODE-AMPLITUDE 
VELOCITY MODE-AMPLITUDE 


A = 
B = 
C = 
FOR 
A = 
B = 
C = 
FOR 
A = 
B = 
C = 


DIMENSION ANR(IO) *DATA(2500) 

COMMON/C OMP1/ QP1*QP2*QP3*QP4*QP5*QP7 
CO MM ON/CO MSI/ Q1 * Q2 * GP1 * GT 1 

COMMON/ FLO DA/ NEQ *UE * AM3* T AU* GAMMA *PI*LIN*TN(10) * I POP 
COMMON/PLTDA/ TARY(lOO) * B ARY ( 10 ► 100 ) 

C0MM0N/C0M2/ B ( 10 ) * C (10 ) * 3R1 ( 10 ) * BR2 ( 10 ) * BR3 ( 10 ) * BS ( 10 * 1 10 > *A(10) 
H * HD2 * HD6* HD8 
T2(l*ln*10)*T3 (4*10*10*10) 

X ( 1 1 ) * IPX(4) 

TSTART.TSTOPfTLYMCY 


400 

402 

403 

405 

406 

407 

408 
410 
430 


C0MM0N/C0M3/ 
C0MM0N/C0M4/ 
C0MM0N/C0M5/ 
C0MM0N/C0M6/ 
FORMAT (8110) 


431 

432 

433 
450 
460 


FORMAT 

FORMAT 

FORMAT 

FORMAT 

FORMAT 

FORMAT 

FORMAT 

FORMAT 

l 

FORMAT 

FORMAT 

FORMAT 

FORMAT 

FORMAT 


INTEGRAL FROM 0 TO 1 OF F(X)'*/ ) 
X*SIN(N*PI*X) *COS(L*PI*X) ' ) 
SIN(N*PI*X>*SIN(M*PI*X)*COS(L*PI*X) ’ ) 
C0S(N*PI*X)*C0S(M*Pl*X)*C0S(L*PI*X) * ) 

X*COS(N*PI*X) *SIN(M*PI*X) *C0S( L*PI*X) * ) 
X*SIN(N*PI*X)*SIN(M*PI*X)*SIN<L*PI*X) * ) 


STEP SIZE = * * F5 . 3 * 


(10X* 'OUTPUT FORMAT 
(10X* 'K=0 IS F ( X ) = 

(10X*'K=1 IS F(X> = 

(10X* 'K=2 IS F(X) = 

(10X* *K=3 IS F ( X ) = 

(10X* 'K=4 IS F ( X ) = 

(BF10.0) 

( 1H1 * 10X* 'SPACE INTEGRALS 

4X * * L =' *12./) 

(/.8X*2HN=* 15* 9110* ) 

( ' K M ' / ) 

(215* 10E10.4) 

(1H1*//*10X* 'DIVERGENT SOLUTION' *//) 

(/. 10X*5HTAU= * F10 . 5 * 5X * 6HNRAR = * FI 0 . 5* 5X * 4HUE = 

15X * 7HGAMMA= . F10 . 5 . 5X * 3 OHP I NI T I AL = *F10.5* 

2//* 10X.6HTIME= *F10.5*5X*6HB(N)= * ElO • 4 * 5X * 6HC < N ) = *E10.4) 

800 FORMAT (5E15.8) 

420 FORMAT (F10. 0*2110) 


* F10 . 5* 


C 

C 

C 


READ SPACE INTEGRALS FROM FILE 2 
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00 200 L=1 * 10 
K = 1 

READ (2*800) ( T2 ( K * N* L) * N=1 * 10 ) 

DO 210 K = l»4 
DO 220 M=1 * 10 

220 READ (2*800) ( T3 ( K * N * M * L > * N=1 ► 10 ) 

210 CONTINUE 
200 CONTINUE 

READ (2*800) HI 

READ INPUT DATA (EXCEPT COM3. PARAMETERS AND INITIAL DISTURBANCE 
FIRST DATA CARD 

NEQ = NO. OF TERMS IN EXPANSIONS 

NX = NO. OF X/L AT WHICH FLOW FIELD CALCULATED 

LIN = 1 TO CALCULATE LINEAR RESULTS 

IPLOT = 1 TO PLOT any output 

INPT = 1 TO WRITE THE SPACE INTEGRALS READ FROM FILE 2 
SECOND DATA CARD 

LC1 = 1 TO WRITE C ( N ) AND B(N) 

LC2 = 1 TO PLOT B ( N ) 

= 4 NO PLOT OF 8 (N) 

LC4 = NUMBER OF TERMS TO BE PLOTTED 
LC5 = INCRIMENTAL INDEX BETWEEN TERMS TO BE PLOTTED 
THIRD DATA CARD 

LP1 = 1 TO CALCULATE U AND P 
LP2 = 1 TO WRITE U AND P 
LP3 = 1 TO PLOT P 

= 4 NO PLOT OF P 

LP4 = NO. OF X/L AT WHICH P OR U TO BE PLOTTED 
FORTH DATA CARD 

NT AU = NO. OF TAU TO BE RUN 
FIFTH DATA CARD 

UE = EXIT MACH NUM8ER 
GAMMA = SPECIFIC HEAT RATIO 
EPS = AMPLITUDE PRECENT error 
SIXTH DATA CARD 

TBEGIN = TIME TO START COMPUTATION OF FLOW VARIABLES AND 
TO START OUTPUT 
TEND = STOP TIME 

TLMCY = START TIME OF LIMIT CYCLE CHECK 
DELTAT TIME DELTA FOR OUTPUT OF LIMIT CYCLE OSCILLATIONS 
TSMPI = START TIME FOR PRESSURE MODE-AMPLITUDE PLOT 
DELPT = TIME DELTA FOR PLOT OF P MODE-AMPLITUDE 
SEVENTH DATA CARD 

X ( I > = AXIAL LOCATION AT WHICH FLOW FIELD IS TO BE CALCULATED 
EIGHTH DATA CARD (USED ONLY IF LP3*4) 

IPX(I) = INDEX OF X ( I ) FOR WHICH PRESSURE IS TO BE PLOTTED 

READ (5*400) NEQ*NX, LIN* IPLOT* INPT 
READ (5*400) LC1*LC2* LC4.LC5 
READ (5*400) LP1*LP2*LP3*LP4 
READ (5*400) NT AU 
READ (5*410) UE*GAMMA*EPS 

READ (5*410) T3EGIN. TEND * TLMCY . DELTAT* TSMPI * DELPT 
READ (5*410) (X(I)*I=1.NX) 

IF (LP3.EQ.4) GO TO 100 
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READ (5*400) ( IPX < I > ♦ 1=1 » LP4) 

100 CONTINUE 

WRITE SPACE INTEGRALS IF INPT = 1 

IF (INPT.NE.l) GO TO 110 
DO 700 L=1 * NEQ 
WRITE (6*430) HI * L 
WRITE (6*402) 

WRITE (6*403) 

WRITE (6*405) 

WRITE (6*406) 

WRITE (6*407) 

WRITE (6*408) 

WRITE (6*431) ( I * I = 1 * NEQ) 

WRITE (6*432) 

M = 0 
K = 1 
J = 0 

WRITE (6*433) J*M* (T2(K*N*L) *N=1*NEQ) 

DO 710 K=l*4 
DO 720 M=1*NEQ 

720 wRITE (6*433) K * M » ( T3 ( K * N * M* L) * N=1 * NEQ) 

710 CONTINUE 
700 CONTINUE 
110 CONTINUE 
C CALL PLOT SUBROUTINE IF IPLOT = 1 
C 

IF (IPLOT. NE.l) GO TO 600 
CALL PLOTS (DATA(l) .2500*3) 

600 CONTINUE 
C 

C CALCULATION OF SOME TERMS USED IN SOLUTION OF ODES 
C 

P2 = 3.14159*3.14159 

UP1=6.28318*UE 

QP2- ( GAMMA-1 • )*UE/2. 

QP 3= (GAMMA-1. ) *P2 
QP4=2.*P2 
GiP5=GAMMA*UE 
QP7=P2 

Q1 = .5* (GAMMA-1. )*UE/GAMMA 
GP1 = GAMMA + 1. 

GT1 = GAMMA*. 5* (GAMMA-1 .) *UE 
Q2 = -Q1*.25*GP1/GAMMA 
C 

DO 1000 KTAU =1*NTAU 
C 

C READ COMBUSTION PARAMETERS 
C 

READ (5*420) TAU*NNB*LCUT 
READ (5*410) (ANR(I) *I=1*NNB) 

C 

LTEMP = TAU/.05 + .01 
H = TAU/LTEMP 
HD2 = H/2. 

HD6 = H/6. 
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HD8 = H/8. 

LTR = (TAU/H) + 1.01 
DO 2000 KK=1*NN3 . 

READ (5*400) NPI 
ANB= ANR(KK) 

DO 510 LLL =1*NEQ 
TN(LLL) = ANB 

IF (LLL.GT.LCUT) TN(LLL) = 0.0 
IF (TAU.LT.l) GO TO 510 
IF (LLL.NE.l) TN(LLL) = 0.0 
510 CONTINUE 

DO 3000 KKK=1*NPI 
READ (5*420) PI*IPOP 
T5TART = TBEGIN 
TSTOP = TEND 
TSMP = TSMPI 
KMT = 0 

KMTS = DELPT/H + 1.01 

KPLT = 2 

K2 = 2 

LGO = 2 

LOUT = 1 

CALL START (LTR*TX*H> 

KONTRL = 2 
L = LTR 

TSTOP1 = TSTOP + .10 
340 IF (TX.GT.TSTOP1) GO TO 130 
IF (L.NE.101) GO TO 140 
LTMP = 102 - LTR 
DO 150 L=1*LTR 
DO 160 I=1*NEQ 
160 BS(I*L> = BS ( I * LTMP ) 

150 LTMP = LTMP + 1 
L = LTR 
140 CONTINUE 

TEST = A8S(TX - TLMCY ) 

IF (TEST. LT. 0.03) K2=l 
IF (K2.NE.1) GO TO 320 
IF (LGO.EQ.l) GO TO 320 
PHIO = 0.0 
DO 900 I=1*NEQ 
900 PHIO = PHIO + 3(1) 

CALL TREND ( TEST * PHI 0 * LGO * EPS > 
IF (LG0.EQ.2) GO TO 370 
TSMP = TX 
TSTART = TX 
TSTOP = TX + DELTAT 
TST0P1 = TSTOP + .10 
370 CONTINUE 
320 CONTINUE 

CHECK = ABS(TX-TSTART) 

IF (CHECK.LT. 0.04) KONTRL = 1 
IF (KONTRL. NE.l) GO TO 330 



IF (LC1.NE.1) GO TO 500 
CALL WOUT1 ( H * T X ) 

500 IF (LC2.EQ.4) GO TO 501 

IF (KMT.GT.KMTS) GO TO 501 
CHK1 = ABS(TX-TSMP) 

IF CCHKl.LE.0.04) KoLT = 1 
IF (KPLT.NE.l) GO TO 501 
KVIT = KMT + 1 
TARY(KMT) = TX 
DO 504 KM=1*10 
BARY(KM*KMT) = 3(KM) 

504 CONTINUE 

IF (KMT.NE.KMT5) GO TO 501 
CALL POUT ( LC4 * LC5 * KMT ) 

KPLT = 2 

501 IF (LP1.NE.1) GO TO 502 

CALL FLOW (NX*H*TX*LP2*LP3*LP4*LP5*L0UT) 

502 CONTINUE 

IF (L0UT.EQ.2) GO TO 3000 
330 CONTINUE 
L = L+l 
TX = TX + H 
L.DO = L-LTR 
LD1 = LDO + 1 
DO 180 I=lrNEQ 
BR1(I) = BS(I*LQ0> 

BR3(I) = BS(l'LDl) 

180 BR2(I) = (BR1<I)+BR3U) )/2. 


CALL RUNG (NEQ) 


300 


DO 300 1=1 • NEQ 
BS ( I * L ) = B(I) 

CHK1 = B ( I > 

CHK2 = C ( I ) , 

IF (CHK1 .LT. 10 . 0 .AND* CHK2.LT. 10.0 ) GO TO 300 

WRITE (6*450) 

WRITE (6*460) TAU * AnB * UE * GAMMA * PI * TX * CHK1 * CHK2 
GO TO 130 
CONTINUE 
GO TO 340 


130 CONTINUE 
3000 CONTINUE 
2000 CONTINUE 
1000 CONTINUE 
STOP 
END 
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400 

410 

420 


SUBROUTINE START ( LTR » TX » H ) 

EXTERNAL POFX 

COMMON/FLODA/ NEQ * Ue . ANB . TAU. GAMMA .PI . LIN. T < 10 ) * IPOP 
COMMON/COV12/ B(10>.C(10>»3R1(10) .RR2(10) »BR3(10) » 8S ( 10. 110 > » A ( 10 ) 
FORMAT (1H1.//.10X. ’INITIAL PRESSURE DISTRIBUTION* . //) 

FORMAT ( 12X» ' Z * . 9X 'P*»/) 

FORMAT ( 5X * 4F10 . 5) 

TX = -TAU 
DO 1 N=1 » NEQ 
A(N) = 0.0 
C ( N ) = 0.0 
B ( N) = 0.0 


DO 100 L =1»LTR 
TX = TX + H 
DO 110 N = 1 * NEQ 
110 8S ( N » L ) = 0.0 
100 CONTINUE 

TX = TX - H 
DO 120 1=1# NEQ 
B ( I ) = POFX <I»PI#IP0P) 
C < I ) = 0.0 
A ( I ) = -B( I) /GAMMA 
120 BS ( I » LTR ) = B ( I ) 

WRITE ( 6 f 400 ) 

WRITE (6.410) 

X = 0.0 

150 SUMB =0.0 


DO 140 1=1. NEQ 
ARG = 3.14159*X*I 
Cl= COS ( ARG ) 

SUMB = SUMB + B ( I ) *C1 
140 CONTINUE 
P = SUM 3 

WRITE (6.420) X.P 
IF (X.GE.1.0) GO TO 200 
X = X + .1 
60 TO 150 
200 CONTINUE 
RETURN 
END 


REAL FUNCTION POFX (I.PI.IPOP) 
IF (IPOP.EQ.il) GO TO 1 

CONTINUOUS WAVE IN IPOP MODE 
POFX = 0.0 

IF (IPOP.EQ.I) POFX = PI 
GO TO 2 

1 continue 

discontinuous il wave 

C = 2 • *PI 
A = 1.5708*1 
POFX = C*SIN(A)/A 

2 CONTINUE 
RETURN 
END 
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SU3R0UTINE TREND ( TEST . PH 1 0 * LSO * EPS) 
DIMENSION PH I MAX ( 6 ) 

IF (TEST. GT. 0.03) GO TO 1 

K = 1 

M=1 

1 IF (M.NE.l) GO TO 10 
IF (PHIO.LE.O) GO TO 4 
PHIBIG = PHIO 

M=2 

GO TO 2 

10 IF (PHIO. LE. PHIBIG) GO TO 3 . 
PHIBIG=PHIO 
60 TO 2 

3 IF (PHIO.GT.O) GO TO 2 
SIGN = PHI0+PHIM1 

IF (SIGN.GT.O) GO To 2 
PHIMAX(K)=PHIBIG 
M = 1 
K = K + 1 

2 PHIM1=PHI0 

IF (K.LE.4) GO TO 4 

AV1=0.0 

AV2=0.0 

DO 5 I=li2 

AV1 = AVI + PHI MAX ( I ) 

IP2= 1+2 

5 AV2 = AV2 +PHIMAX ( IP2) 

K - 1 

DELTA = A3S( ( AV2-AV 1 ) /2 . 0 ) 

CHECK= EPS*AV2/2.0 
IF (DELTA. GT. CHECK) GO TO 4 
L60=l 
GO TO 6 

4 LG0=2 

6 CONTINUE 
RETURN 
END 



SUBROUTINE FLOW ( NX . H » TX . LP2 . LP3 » LP4 . LP » LOUT ) 

COMMON/FLOD'A/ NEQ » UE » ANB » TAU. GAMMA » P I • LIN * T ( 1 0 ) - 

COMMON/ C0M2/ B ( 10 ) * C < 10) » BR1 ( 10) » r'-.R2( 10) » 3R3( 10) * BS( 10» 110) » A ( 10 » 
COMMON/COM5/ X(ll) . IPX (4) 

C0MM0N/C0M6/ TSTART.TSTOP.TLYMCY 

COMMON/C 0M8/ A3C(303> » ORDP ( 1 1 » 303 > » ORDU ( 1 1 » 303 ) 

400 FORMAT ( 1 H 1 # /.lOX.’FLOW PRARAMETERS* »10X.3HZ= »F6.3»/ ) 

410 FORMAT < 10X» 'LINEAR= • » 1 2 * 9X . « NEG = *.I2» 9X.'N= **F7.5. 9X»»TAU= * 

1 »F7.5./» 10X. 'EXIT MACH = ' » F5 . 3 » 3X . ’ GAMMA =*»F5.3. 4X»*PINITIAL =' 

2 »F5.4 ./) . . 

420 FORMAT (11X» 4HTIME.3X. BHPRESSURE » 3X » 8HVEL0C ITY »/ ) 

430 FORMAT ( 10X » F7. 3. 1 1E10 . 4 ) 

TEST = ABS ( TX-TST ART ) 

IF (TEST. GT. 0.03) GO TO 1 
K = 1 

1 CONTINUE 

DO 110 N=1*NX 
A1 = 3. 14159*X ( N ) 

VEL=X(N) *UE 
SUMA = 0.0 
SUMS = 0.0 
SUMC = 0.0 
SUMU = 0.0 
00 120 1=1# NEQ 
TA= A1*I 
ST = SIN(TA) 

CS = COS(TA) 

SUMA = SUMA + A ( I ) *CS 
SUMS = SUMB + B ( I ) *CS 
SUMC = SUMC + C ( I ) *ST 
SUMU = SUMU + C(I)*I*3.14159*ST 
120 CONTINUE 

CALL PRMTRS ( N . K » SUMA t SUMB » SUMC » SUMU. VEL ) 

110 CONTINUE 

A3C ( K ) = TX 

IF (TX.LT.TSTOP) SO TO 300 
LOUT = 2 
KSTOP = K 

IF (LP2.NE.1) GO TO 200 

DO 310 J = 1 » NX 

KOUNT = 44 

DO 220 L=1»KST0P 

IF (K0UNT.NE.44) GO TO 210 

WRITE (6.400) X(J) 

WRITE (6.410) LlN»NEQ»ANB.TAU*UE .GAMMA. PI 
WRITE (6.420) 

KOUNT = 1 
210 CONTINUE 

WRITE (6.430) ABC (L ). ORDP ( J * L ). ORDU ( J » L) 

KOUNT = KOUNT + 1 
220 CONTINUE 
310 CONTINUE 

200 IF (LP3.EQ.4) GO TO 300 

CALL P0UT2 (LP3.LP4. KSTOP. NX) 

300 CONTINUE 
K = K + l 
RETURN 
END 
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SUBROUTINE P0UT2 < LP3 * LP4 * KSTOP » NX > 

COMMON/FLODA/ NEQ * UF > ANB * T AU* GAMMA * PI » LI N * T ( 10 ) 

COMMON/ CO VI 3/ H 
C0MM0N/C0M5/ X(11)*IPX(4) 

C0MM0N/C0M8/ A9C(303) * ORDP ( 11 * 303 ) * ORDU( 11 * 303) 

C0MM0N/C0M9/ 0RQ(303) 

CALL PLOT (0.0*2.0*-3> 

CALL PLOT (0.0*11.0*3) 

CALL PLOT (1.0»0.5»-3) 

TERMS = NEQ 

NPT = KSTOP 

J1 = NPT + 1 

J2 = NPT + 2 

SIZE = 0.10*NPT 

CALL SCALE (A3C*SIZE*NPT*1) 

DO 1 J=1*LP4 
DO 4 I=1*NX 
ICHK = IPX(J) 

IF (ICHK.NE.I) GO TO 4 
Z = X 1 1 ) 

□0 100 M = 1 * KSTOP 
100 ORD(M) = ORDP ( I * M ) 

GO TO 110 
4 CONTINUE 
110 CONTINUE 

CALL SCALE ( ORD *4.0* NPT » 1 ) 

IF (J.NE.3) GO TO 2 
DELX = SIZE + 4.0 
CALL PLOT ( DELX * -6 • 3 * ~3 ) 

2 IF ( J.E3.2.0R. J.EQ.4) GO TO 3 

CALL SYMBOL ( 2 . 90 * 1 . 80 * 0 . 10 * 32HN0RMALI ZED PRESSURE "TIMF HISTORY* 
1 0 • U * 32 ) 

CALL SYM30L (2.30*1.55*0.10* 3HN =*0.0*3) 

CALL SYMBOL (3.30*1.55*0.10* 4HTAU=* 0.0*4) 

CALL SYMBOL (4.50*1.55*0.10* 4HUE =*0.0*4) 

CALL SYMBOL (5.60*1.55*0.10* 6HGAMMA=* 0.0*6) 

CALL SYMBOL (2.30*1.30*0.10* 5HNES =*0.0*5) 

CALL SYMBOL (3.30*1.30*0.10* 3HH= *0.0*3) 

CALL SYM30L (4.50*1.30*0.10* 4HPI =*0.0*4) 

CALL NUM3ER ( 2 . oO * 1 . 55 * 0 . 1 0 * ANB * 0 . 0 * 4 ) 

CALL NUM3ER ( 3 . 75 * 1 . 55 * 0 . 1 0 * T AU» 0 . 0 * 4 ) 

CALL NUMBER (4.95*1.55*0.10* UE*0.0*3) 

CALL NUMBER ( 6 . 25 * 1 . 55 * 0 . 10 * GAMMA * 0 . 0 * 3 ) 

CALL NUMBER ( 2 • 80 * 1 • 30 * 0 • 1 0 * TERMS * 0 • 0 * -1 ) 

CALL NUMBER ( 3. 55 * 1 . 30 * 0 . 1 0 * H * 0 . 0 * 3 ) 

CALL NUM3ER ( 4 . 95 * 1 . 30 * 0 . 1 0 * Pi * 0 . 0 * 3 ) 

3 IF ( J.EQ.2.0R. J.EQ.4) DEL Y = 5.3 
IF ( J.EQ. l.OR. J.EQ.3) DELY = 4.0 
DELX = 0.0 

IF (J.EQ. l.OR. J.EQ. 3) DELX = 2.0 
CALL PLOT ( DELX * DELY * ”3 ) 

CALL SYM30L ( 1 . 80 * -1 . 70 * 0 . 1 4 . 4HX/L=# 0.0 *4) 

CALL NUMBER ( 2. 40 * -1 . 70 * 0 . 14 * Z*0.0'3) 

CALL FACTOR (0.788) 

CALL AXIS (0. 0*0.0* 4HTI ME * -4 * SIZE * 0 . 0 * ABC ( J1 ) t ABC ( J2 ) ) 
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CALL AXIS (0.0.-2.0, 8HPRESSURE , 8 , 4 . » 90 . , ORD < J1 ) , ORD ( J2) ) 
CALL PLOT ( 0 . 0 , -2 • 0 » -3 ) 

CALL LINE (ABC, ORD, NPT, 1,1,1) 

CALL FACTOR (1.0) 

1 CONTINUE 
RETURN 
END 


SUBROUTINE POUT (NMDE» ISP.NPT) 

DIMENSION COEF(IOO) 

COMMON/FLODA/ NEQ » UE > ANB » TAU * GAMMA . PI , LIN» T ( 1 0 ) 

COMMON/PLTDA/ TIM(lOO), BS(10»100) 

C0MM0N/C0M3/ H 
EON = NEO 
J1=NPT+1 
J2=NPT+2 

CALL SCALE ( TIM, 4 • 0 » NPT » 1 ) 

CALL PLOT (0.0, 2. 0,-3) 

CALL PLOT (0. 0.11.0.3) 

CALL PLOT (0.0. 0.5. -3) 

KOUNT = 1 

DO 110 1=1 . NMDE. ISP 
L = 1 

DO 120 K = 1 » NPT 
COEF ( K ) =BS ( I » L ) 

120 L=L+1 

CALL SCALE ( COEF . 2 . 0 . NPT . 1 ) 

IF (KOUNT. NE.l) GO TO 1 

CALL SYMBOL (3. 0.1. 3»0,10. 27hT I ME DEPENDENT COEFFICIENTS . 0 . 0 . 27) 
CALL SYMBOL (3. 5.1. 1.0.10. 14HTERM EXPANSI ON » 0 . 0 . 14 ) 

CALL SYMBOL ( 3 . 0 . 0 . 9 . 0 . 10 . 29HTAU= NBAR= H=.0.0»29) 

CALL SYMBOL (3. 0.0. 7.0. 10* 20HGAMMA= MACH= ,0.0.20) 

CALL NUMBER (3.0,1.1,0.10, EQN.O.C.-l) 

CALL NUMBER ( 3 . 5» 0 . 9 » 0 . 10 , TAU, 0 . 0 » 3 ) 

CALL NUMBER (4.9,0.9,0.10»ANB»0.0,3) 

CALL NUMBER ( 6. 0 » 0 . 9 » 0 . 10 , H » 0 . 0 , 3 ) 

CALL NUMBER ( 5. 1 . 0 . 7 » 0 . 1 0 , UE ,0.0,3) 

CALL NUMBER (3.7.0.7,0.10, GAMMA ,0.0,3) 

CALL PLOT (2. 0,0. 5, -3) 

1 Y = 2»5 

CALL PLOT (0.0.Y.-3) 

CALL AXIS ( .0, .0,4HTIME,-4»4.,0.0,TIM(J1) ,TIIM(J2) ) 

CALL AXIS (0.0, -1.0, 5HBN( T ) ,5» 2.0, 90.0 , COEF ( J1 ) » COEF ( J2 ) ) 

CALL SYMBOL ( 4 . 5, 0 . 0 , 0 . 10 , 2HN=» 0 . 0 ,2) 

TERM = I 

CALL NUMBER ( 4 . 8 » 0 . 0 » 0 . 10 » TERM » 0.0,-l> 

CALL PLOT (0. 0,-1. 0.-3) 

CALL LINE (TIM.COEF. MPT, 1,0,1) 

CALL PLOT (0.0, 1.0, -3) 

IF (KOUNT. NE. 3) GO TO 20 
CALL PLOT ( 8. 0,-8. 0,-3) 

KOUNT =1 
GO TO 110 

20 KOUNT = KOUNT + 1 
110 CONTINUE 

CALL PLOT (8. 0,0. 0,-3) 

CALL PLOT (0.0,0.0,999) 

RETURN 

END 
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Section of the Program. Used to Solve the Nonlinear Wave Equation 


SUBROUTINE RUNG (NEQ) 

INTEGRATION OF SECOND ORDER ODE WITH RETARDED VARIABLE. 
USE WITH SECOND ORDER WAVE EQUATION ANALYSIS. 


100 

110 

120 

130 

140 


EXTERNAL EQTN 

DIMENSION R ( 10 > 4) *BB( 10 ) * BPB ( 10> * RBI ( 10) * BPB1 ( 10) , 

C0MM0N/C0M2/ BP (10) > B ( 1 0 ) *Rl < 10 ) *R2« 10 ) *R3 ( 10 ) * BS< 10 . 1 10 > * DM ( 10 ) 

COMMON/COM3/ H*H2*H06*H8 

DO 100 1 = 1 * NEQ 

R(I»1> = H*EQTN(I. R 1 * 3 * BP ) 

BPB(I) = BP(I>+ R(I.l)/2. 

BB(I) = B ( I ) + H2*BP ( I ) + H8*R(I*1) 

DO 110 1 = 1* NEQ 

R(I*2) = H*EQTNtI. R2 *3B*BPB) 

BPB1 ( I ) = BP ( I ) + R ( I * 2 ) /2 . 

BB1(I)= B<I) + H2*3P ( I ) +H8*R (1*1) 


DO 120 1 = 1 * NEQ 

R ( I • 3) = H*EQTN(I. R2 .BB1.BPB1) 

BPB(I)= BP ( I ) + R ( I * 3 ) 

BB(I) = B(I) + H*3P ( I ) ♦H2*R(I*3) 

DO 130 1 = 1 * NEQ 

R ( I » 4 ) = H*EQTN ( I * R3 *BB*BPB> 

DO 140 1=1. NEQ 

B(I) = H*(BP(I)+(R(I*1)+R(I*2>+R(I*3) )/6.) + Bd) 
BP(I)= (R( I * 1)+2.*(R(I *2) *Rd *3> )+R(I *4) )/6. +BP(I) 
RETURN 
END 


C 

C 

c 


REAL FUNCTION EQTN ( L * YPR * Y * YP > 
SECOND ORDER WAVE EQUATION 


110 

1 

100 


DIMENSION Y(10) *YP(10> *YPR(10> 

COMMON/COMP1/ Q1*Q2*Q3*Q4«Q5*Q7 
COMMON/COM4/ T2(1*10*10)*T3(4*10*10*10) 

COMMON/FLOO A/ NEQ * UE * AN8* TAU* GAMMA » PI * LIN* T <10 ) 

Di = -L*L*Q7*Y(|_) - Q5*(YP(L) - T(L>*(YP(L> - YPR(L) ) ) 


SUM=.0 

DO 100 N=1 * NEQ 

51 = Q1*N *YP(N)*T2(1*N*L) 

52 =((-l)**(N+L ) ) * YP < N ) *Q2 
SUMi=. 0 

IF (LIN.EQ.1) GO TO 1 
DO 110 M=1 * NEQ 

53 = Q3* M*M *YP(N)*Y(M)»T3(2*N*M*L) 

54 = Q4* N*M *Y(N)*YP(M)*T3( 1*N*M*L) 
SUMl= SUM1+S3-S4 

CONTINUE 

SUM = SUM + SUM1 ♦Sl-SB 
EQTN = 01 +2.*SUM 
RETURN 
END 
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SUBROUTINE PRMERS (N,K,SUMA,SUMB ,SU1C,SUMU,VEL) 

SUBROUTINE FOR CALCULATING FLOW PARAMETERS FOR WAVE EQUATION 

COMMON/FLODA/ NEQ , UE , ANB , TAU , GAMMA ,PI , LIN , T (ID) 

COMMON/ COM5 / X(ll) ,IPX(4) 

COMMON/COM8/ ABC(303) ,0RDP(11,303) 2 0RDU(11j3O3) 

ORDU(N,K) = -SUMU 
IF (LIN.EQ.l) GO TO 1 

ORDP(N,K) = GAMMA*(SUMB*(SUMB-2.)-+6UMU*(2.*VEL-SUMU))/2. 

GO TO 2 

1 ORDP(N,K) = GAMMA* (-SUMB + VEL*SUMU) 

2 CONTINUE 
RETURN 
END 


SUBROUTINE WOUTl (H.TX) 

COMMON/FLODA/ NEQ * UF * ANIB * TAU# GAMMA * PI * LI N * T ( 10 ) 

C0MM0N/C0M2/ 8 < 10 ) . C ( 10 ) . 3R 1 ( 10 ) . RR2 ( 10 > * 9R3 ( 10 ) . 3S < 1 0 . 11 0 > . A ( 10 > 
C0MM0N/C0M6/ TST ART , TSTOP * TL YMC Y 
420 FORMAT < 3X . F7 . 3 , 10E10 . 4 ) 

430 FORMAT (1H ) 

440 FORMAT ( 1H1 * 1 OX * • TI ME DEPENDENT COEFFICENTS OF THE 

1 » * NONLINEAR WAVE EQUATION PHI = A(T)*COS(n*PI*Z> • */ 

( 5X * ’TIME 
AP7 

(5X» ’TIME 
A7 


451 FORMAT 
1 AP6 

452 FORMAT 

1 A6 
410 FORMAT 

1 

2 


API 

AP2 


AP3 


APB 


AP9 


AP10 * ) 

A1 

A2 


A3 


AB 


A9 


A10’ ) 


AP4 


A4 


AP5 


A5 


( 10X * ’ LINEAR= ' * I2*9X* ’NEQ= **I2* 9X*’N= »*F7.5* 9X* 


* F7. 5* / * 10X ► 'EXIT 
*F5.4./) 


TAU= 


MACH= ' * F5 * 3* 3X * ’GAMMA ='*F5.3* 4X*’PINITIAL = » 


TEST = ABS(TX-TSTART) 

IF (TEST. GT . 0.030) GO TO 10 
K = 16 

10 IF (K.NE.16) GO TO 2 
WRITE ( 6 i 440 ) 

WRITE (6*410) LIN.NEQ* ANB*TAU*UE*GAMMA*PI 
WRITE ( 6 * 452 ) 

WRITE (6*451) 

K = 1 

2 WRITE (6*430) 

WRITE (6*420) TX* (C ( I ) * 1=1 * NEQ) 

WRITE (6.420) TX» <B ( I ) * 1 = 1 . NEQ) 

K=K + 1 

RETURN 

END 
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Section of the Program Used in the Analysi s of Large Arrylitude 
Oscillations 


SU3R0UT I N£ RUNG <NEQ) 

integration of first order ode with retarded variable 

USE WITH LARGE AMPLITUDE ANALYSIS 

TD^cJi“ , Bui!"<10..=.RH10..BR 2 (10..B R 3<10 1 .BS<10.nO».A 1 10. 

COMMON/COM3/ H»HD2»HD6 

CALL EQTN ( A t B * C # BR 1 • R A 1 # RBI » RC 1 ) 

DO 100 I =1*NEQ 

A1(I) = A(I) + HD2*RAld) 

Bid) = 8 d ) + HD2*RB1(I) 

C1(I> = CCI) + HD2*RC1(I) 

100 CALL^EQTN <A1»B1*C1»BR2*RA2»RB2»RC2) 

DO 110 I=1»NEQ 

A2CI) = A(I) + HD2*RA2(I) 

B2CI) = BCD + HD2*RB2 C I ) 

C2CI) = CCD .+ HD2 + RC2CI) 

110 CAL^EQTN (A2*B2*C2»BR2*RA3*RB3»RC3) 

DO 120 I=1>NEQ 
A 1 ( I ) = ACD + H*RA3d) 

Bl(I) = BCD + H*RB3 C I ) 

Cl c I ) = CCI) + H*RC3d> 

120 CALL. 1 EQTN C A 1 » B1 * Cl . BR3 * RA4 1 R34 » RC4 ) 

MD 3 = JfD N + Q H06dRAl(l)>RA4Cl) + 2.MRA2CIDRA3CI))) 

BCI) = BCD + HD6*CrB1<I)+RB4CI)+2.*CRB2 I +RB3I) 

CCI) = CCD + HD6* C RC1 Cl ) +RC4C I) +2 . *C RC2 Cl ) RC 
130 CONTINUE 
RETURN 
END 
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SUBROUTINE PRMTRS (N,K,SUMA,SUMB,SUMC,SUMU/VEIj) 

SUBROUTINE FOR CALCULATING FLOW PARAMETERS FOR LARGE AMPLITUDE WAVES 

COMMON/FLOM/ NEQ,UE, ANB,TAU, GAMMA, PI s LIN,T(LQ'3) 

COMMON/COMSl/ Q1,Q2,GP1,GT1 
COMMON/COM5/ X(ll) ,IPX(4) 
common/com8/ abc( 303 ) jOrdpCh^os) ,ordu(h, 303 ) 
ordp(n,k) = sumb 

IF (LIN.EQ.l) GO TO 1 

ORDU(N,K) = SUMO + (Q1 + Q2*SUMB) *SUMB*X(N) 

GO TO 2 

1 ORDU(N,K) = SUMC + Q1*SUMB*X(N) 

2 CONTINUE 
RETURN 
END 


420 

430 

440 


SUBROUT I ME W0UT1 (H,TX) 

COMMON/FLODA/ NEQ . UE * ANB * TAU. GAMMA , PI * LIM * T ( 10 > 

C0.MM0N/C0M2/ B<10) *C(l0) * BR1 < 10 ) ► ?,R2 < 1 0 ) * BR3 ( 10 > *BS(10*110)*A(10> 
COMMON/COMSl/ 31 * 02 * GP1 * GT1 
COMMON/ C0M6/ TSTART,TSTOP*TLYMCY 
FORMAT (3X*F7.3*10E10.4) 

FORMAT (1H ) 

FORMAT (lHlflOX* 'TIME DEPENDENT COEFFICENTS OF THE '* 

•LARGE AMPLITUDE SOLUTIONS' */) 


450 FORMAT (5X*'TI.ME 

1 B6 B7 

451 FORMAT (5X»'TIME 

1 C6 C7 

452 FORMAT (5X* 'TIME 

1 A6 A 7 

410 FORMAT < 10X * • LI NEAR= ' * I 2 * 9X * » NEQ = »*I2* 9X*'N= **F7.5* 9X*'TAU= 

1 * F7.5* / • 10X * ' EXI T MACH= ' ► F5 • 3 * 3X * ' GAMMA ='*F5.3* 4X*'PINITIAL : 

2 * F5.4* /) 


B1 


B2 


B3 



B8 


B9 


BIO' ) 

Cl 


C2 


C3 



C8 


C9 


CIO') 

A1 


A2 


A3 



AB 


A9 


A10' ) 


B4 


C4 


A4 


B5 

C5 

A5 


TEST = A3S ( TX-TST ART) 

IF (TEST. GT. 0.030) GO TO 10 
K = 12 

10. IF (K.NE.12) GO TO 2 
WRITE (6*440) 

WRITE ( 6 r 4 1 0 ) LIN*NEO* ANB * TAU* UE * GAMMA * PI 
WRITE ( 6 ► 452 ) 

WRITE ( 6 * 450 ) 

WRITE (6*451) 


K = 1 

2 WRITE (6*430) 

WRITE (6*420) TX * ( A ( I ) * I = 1 * NEQ ) 

WRITE (6*420) TX * ( B ( I ) * 1 = 1 * NEQ ) 

WRITE (6*420) TX* (C(I) *1=] *NE3) 

K=K + 1 

RETURN 

END 
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SUBROUT I ME EQTN ( A * 3 • C * 8R • RA » RB » RC ) 


- T<L>*UE*(B(L> - BR(L> > 


LARGE AMPLITUDE EQUATION 

DIMENSION A(10)»B(10)»C(10)iBR(10).RA(10).RB(10)»RC(10).F1(10)» 

1F2(10) »F3( 10) 

COMMON/COMS1/ Q1 » 02 » GP1 * GT1 

COMMON/FLODA/ NEQ»UE» ANB* TAU* GAMMA »PI »LIN* T ( 10 > 

COMMON/COM4/ T2(l.l0»10)»T3(4»10»10»10) 

DO 100 L=1*NEQ 
PIL = L*3* 14159 
SOI =PIL*C<L> + Q1*R(L> 

50 = -UE*A(L) + SOI 
RO = “GAMMA* ( S01+ UE*9(L)) 

UO = PIL*B (L) /GAMMA - UE*C(L> 

SUMN1 = 0.0 
SUMN2 = 0.0 
SUMN3 = 0.0 
DO 110 N=1*NEQ 
PIN = N*3. 14159 

51 = PIN*T2 ( 1 • N*L) *A (N) 

52 = PIN*T2 (-1 *N*L) *R(N) 

R1 = PIN*T2( 1 r 4 * L 1 * B < N ) 

Ul = PIN+T2(1*L*N)*C(N) 

SUMM1 = 0.0 
SUMM2 = 0.0 
SUMM3 = 0.0 

IF (LIN.EQ.l) GO TO 200 
DO 130 M=1 * NEQ 
PIM = M*3. 14159 

53 = P1M*T3(1»N»M*L)*C(N)*A(M) 

54 = PIM*T3(2»N»M»L)*A(N)*C<M) 

55 = T3(2»N»M»L)*A(N) *A<M) 

56 = PIM*T3(3»N»M»L)*<3(N)*A(M) 

58 = T3(2*N»M»L) *A(n) *3(M) 

59 = T3(2»NfM»L)*(B(N) -BR(N) )*A(M) 

510 = T3(2»N»M#L)*B(N)*B(M) 

511 = PIM*T3(3»N»M»L)*9(N)*B(M) 

512 = S10 - 2.*S11 
R2 = PIM*T3(1»N»M»L)*C(N)*3(M) 

R3 = PIM*T3(2»N.M»L)*3(N)*C(M) 

R4 = T3(1»N*M»L)*C(n)*C(M) 

U2 = PIM*T3 ( 1 » L« M»N) *A (N) *B < M> 

U3 = T3(1»N»L»M)*C(N)*CCM)*PIM 
U4 = PIM*T3( 3* N»L* M) *3 ( N) *C ( M) 

U5 = PIN*T3(4»N»L*M)*3(N)*C(M) 

U6 = T3(l*L»MrN)*B(N)*C(M) 

SUMM1 = SUMM1 + S3 * S4 - UE*S5 
1 + Q2*S12 

SUMM2 = SUMM2 + R2 - GAMMA*R3 “ 

1 - Q2*GAMMA*S12 

SUMM3 = SUMM3 + U2/GAMMA - U3 - 
130 CONTINUE 
200 CONTINUE 

SUMN1 = SUMN1 + SUMM1 + UE*S1 - 


- B<M>*A(N) ) 


+ Ql*(S6+S8> -2.*T(N)*UE*S9 
GT1*R4 - Ql* (GAMMA+S10 - GP1*S11) 


Q1*(U4 “ U5 + U6) 


Ql*S2 
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SUMN2 = SUMN2 + SUMV2 + UE*R1 + &AMMA*Q1*S2 
SUMN3 = SUMN3 + SUMM3 - UE*Ul 
110 CONTINUE 

F1(L) = SO + 2.*SUMN1 
F2tL) = RO F 2 < *SUMN2 
F3(L) = UO + 2 • *SUMN3 
100 CONTINUE 

DO 300 L=1»NEQ 
U8N =0.0 
DO 310 N = 1 » NEQ 
U8 = T2(1»L»N>*F2(N) 

U8M = 0.0 

IF (LIN. EQ. 1 ) GO TO 320 

DO 330 M=1 t NEQ 

U9 = T3(3*N»L»N)*a(N)*F2(y) 

330 U8M = U8V1 + U9 
320 CONTINUE 

U8N = U8N - Q1*U8 - Q2*U8^*2. 

310 CONTINUE 

RA(L) = F1(L) 

R3(L) = F2(L) 

RC <L) = F3(L> + 2 . *IJ8N 
300 CONTINUE 
RETURN 
END 
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